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ABSTRACT 


Is  this  thesis,  we  examine  the  ssefolness  of  the  two-parameter  Born 
inversion  method  of  Clayton  and  Stolt  (1981)  nnder  idealized  conditions.  We 
implement  a  constant  background  procednre  to  reconatroct  line  source 
synthetic  data  which  simnlate  a  stratified  aconstic  medium.  Specifically, 
this  investigation  is  an  extention  of  the  sing*,  interface  work  of  Weglein, 
Violette  and  Eeho  (1985)  for  doable  interface  models.  Since  we  nse  offset 
data  to  recover  one-dimensional  bnlk  modulus  and  density  variations,  this 
procedure  is  termed  a  1.5  dimensional  inversion.  Implicit  in  these 
reconstructed  acoustic  profiles  are  errors  related  to  the  approximate 
inversion.  Phase  shifts  inherent  in  the  modeled  wave  field  cause  these 
inversion  errors  not  only  to  propagate,  but  also  to  increase  with  depth. 
Since  the  model  data  are  exact,  we  are  able  to  examine  the  analytic 
expressions  of  these  inversion  artifacts.  From  this  analysis,  we  develop  an 
implementation  which  minimizes  inversion  errors.  This  it  achieved  by  both 
preconditioning  the  offset  data,  and  judiciously  choosing  a  reference 
velocity  in  accordance  with  inversion  goals.  Both  of  these  procedures 
however  provide  limited  improvement,  and  we  conclude  that  this  constant 
background,  frequency  domain,  method  should  be  enhanced  by  more 
sophisticated  algorithms.  Specifically,  we  suggest  the  implementation  of 
either  a  variable  background  procedure,  or  the  refinement  algorithm  of  Bagin 
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and  Cohen  (1984). 
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PREFACE 


In  seismic  exploration,  an  energy  source  located  at  the  Earth's  surface 
initiates  seismic  waves  which  propagate  through  the  subsurface.  A  portion 
of  this  seismic  energy  is  reflected  by  inhomogeneities  within  the  probed 
medium,  and  this  reflected  response  is  systematically  measured  along  the 
Earth's  surface.  In  general,  the  object  of  seismic  exploration  is  to 
determine  the  nature  of  the  subsurface  from  these  geophysical  data.  That 
is,  given  knowledge  about  both  the  source  and  receivers,  we  hope  to  extract 
subsurface  information.  This  information  depends  upon  both  the  mathematical 
model  used  to  describe  the  medium,  and  the  method  of  inversion.  In  this 
investigation,  we  implement  a  method,  founded  in  inverse  scattering,  for 
recovering  both  acoustic  parameters. 

In  our  method,  we  assume  that  subsurface  parameter  changes  can  be 
modeled  as  deviations  about  a  known  reference  medium.  Thus,  as  the  incident 
wave  field  propagates  through  the  medium,  any  variation  of  the  wave  field 
from  the  reference  wave  field  is  attributed  to  scattering.  Expressing 
reflected  energy  as  these  scattering  interactions  establishes  the 
mathematical  framework  for  recovering  subsurface  information.  Basic  to  any 
mathematical  expression  for  wave  propagation  is  a  physical  model  describing 
the  assumed  nature  of  the  medium.  Typically,  these  expressions  are  partial 
differential  equations  which  include  spatially  variant  coefficients. 
Determining  these  coefficients  determines  the  model  parameters.  As 
mentioned  above,  we  assume  an  acoustic  model  and  recover  the  corresponding 
bnlk  modnlns  and  density  variations. 

Scattering  theory  is  characterized  by  both  the  forward  and  inverse 
problems.  In  the  forward  problem,  we  specify  the  model  parameters  and 


calculate  the  corresponding  wave  response.  This  procedure  has  been  used  in 
an  iterative  scheme  to  recover  subsurface  information/  given  a  data  set,  we 
would:  (1)  calculate  a  forward  response,  (2)  compare  this  synthetic  response 
with  the  observed  response,  and  (3)  minimize  the  resulting  model  parameter 
differences  through  an  iterative  procedure.  Inverse  scattering  however 
seeks  to  directly  determine  the  model  parameters  from  the  observed  data. 
These  methods  are  constrained  by  the  specific  assumptions  which  characterize 
the  model.  In  this  thesis,  we  consider  an  approximate  inversion  scheme  for 
directly  recovering  one-dimensional  acoustic  parameter  variations. 

Since  this  inversion  scheme  is  approximate,  subsurface  phenomena  are 
either  incorrectly  modeled,  or  not  modeled  at  all.  Specifically,  our  method 
is  a  linear  procedure  which  ignores  both  multiple  reflections  and 
transmission  losses.  Furthermore,  all  reconstructed  parameter  changes  are 
contaminated  by  linearization  errors.  We  therefore  focus  our  attention  on 
the  reliability  of  this  direct  inversion  method. 
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INTRODUCTION 


The  goal  of  seismic  inversion  is  to  recover  snbsnrface  properties  from 
surface  reflection  data.  For  a  particular  method,  these  recovered 
properties  depend  on  the  assumed  nature  of  the  probed  medium.  In  this 
thesis,  we  examine  the  nsefulness  of  a  well  known  linear  inversion  procedure 
which  extracts  acoustic  parameters  from  a  stratified  earth.  Before 
discussing  this  method  however,  we  outline  the  development  of  other  related 
techniques . 

Prior  to  the  development  of  direct  inversion  procedures,  structural 
inversion  was  accomplished  by  wave  equation  migration  methods.  These 
methods  have  been  refined  and  are  widely  used  today.  Although  wave  equation 
migration  methods  can  accurately  locate  subsurface  reflectors,  they  can  only 
qualitatively  predict  the  sizes  of  corresponding  impedance  contrasts — hence 
the  term  structural  inversion.  Briefly,  all  wave  equation  migration  methods 
rely  on  specific  imaging  assumptions  to  map  subsurface  reflectivity  by 
backpropagating  the  recorded  wave  field.  Recent  survey  books  discuss  both 
the  theoretical  background  (Morgan,  (1983)),  and  the  available  techniques 
(Berkhout,  (1982)  and  Robinson,  (1983))  of  wave  equation  migration. 

In  a  1979  paper,  Cohen  and  Bleistein  present  and  implement  a  procedure 
for  the  direct  inversion  of  the  acoustic  wave  equation  for  a  two-dimensional 
velocity  profile.  Expressing  the  forward  equation  as  Green's  theorem,  they 
linearize  the  resulting  integral  equation  and  adopt  a  perturbative  approach 
before  inverting.  This  linearization  admits  two  assumptions  inherent  in 
wave  equation  migration:  (1)  all  reflection  events  are  treated  as  primary 
reflections,  and  (2)  a  constant  velocity  medium  lies  above  each  reflector. 
As  demonstrated  by  Cohen  and  Bleistein,  the  linearization  does  not  seriously 
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impair  their  method.  Io  fact,  this  direct  inversion  offers  an  improvement 
over  the  imaging  techniques  of  wave  equation  migration.  In  direct 
inversion,  not  only  are  the  reflectors  imaged,  bat  also  the  sizes  of 
impedance  changes  are  recovered  within  the  linearization  assumption.  In 
effect,  this  inversion  is  a  direct  mapping  procedure  which  assigns  a  value 
of  acoustic  velocity  to  a  specific  subsurface  location.  In  subsequent  work, 
Bleistein  and  Cohen  (1982)  present  a  three-dimensional  inversion  algorithm 
as  well  as  implement  their  method  for  field  data.  Also,  Hagin  and  Cohen 
(1984)  refine  the  previous  method  to  compensate  for  the  effects  of  the 
linearization. 

Both  wave  equation  migration  and  the  approximate  inversion  of  Bleistein 
and  Cohen  utilize  the  conventional  common-midpoint  stacked  section  to 
simnlate  a  coincident  source-receiver  configuration.  The  processing  of  this 
data  however  eliminates  all  information  concerning  the  variation  of  the 
reflection  coefficient  with  incident  angle.  Thus,  these  techniques  cannot 
separate  density  variations  from  velocity  or  bulk  modulus  variations.  As 
shown  by  Clayton  and  Stolt  (1981)  and  Raz  (1981),  the  above  method  of  linear 
inversion  provides  the  mathematical  basis  for  recovering  both  velocity  and 
density  variations.  Both  papers  present  alternate  approximate  inversion 
algorithms  for  separating  the  components  of  impedance.  Since  both 
algorithms  are  formulated  from  Born  theory,  their  underlying  assumptions  are 
identical  to  those  of  Cohen  and  Bleistein.  In  fact,  both  methods  reduce  to 
Cohen  and  Bleistein's  result  for  zero  offset  and  constant  density.  Hanson 
(1984)  has  demonstrated  this  in  his  discussion  of  the  three  methods. 

Although  the  Born  inverse  methods  of  Clayton  and  Stolt,  and  Raz  both 
rely  on  the  angular  dependence  of  the  reflection  coefficient  to  separate 
velocity  and  density  variations,  their  approaches  are  different.  The  method 


of  Clayton  and  Stolt  is  similar  to  that  of  Cohen  and  Bleistein  in  that  the 
data  mast  be  Fourier  transformed  over  both  time  and  the  spatial  coordinates. 
In  contrast,  Raz's  procedure  operates  in  the  time  domain.  In  implementing 
Raz's  method,  Hanson  (1984)  demonstrates  that  this  time  domain  algorithm  has 
greater  stability  than  the  frequency  domain  methods.  According  to  Hanson, 
this  improved  atability  derives  from  the  operator's  ability  to  accurately 
record  amplitude  variations  with  offset.  (Frequency  domain  methods 
manipulate  the  entire  wave  field  at  once,  and  thus  an  accurate  reference 
velocity  is  necessary  for  stability.)  Hanson  admits  however  that  the 
careful  picking  of  events,  critical  to  Raz's  method,  is  a  tedious  chore. 
Thus,  although  Raz's  method  leads  to  more  accurate  results,  there  is  of  yet 
no  automated  method  of  implementation. 

In  this  investigation,  we  reformulate  and  implement  the  frequency 
domain  method  of  Clayton  and  Stolt  for  synthetic  data  representing  a 
stratified  acoustic  medium.  Since  we  are  reconstructing  one-dimensional 
parameter  variations,  we  require  only  a  single  common  shot  gather.  To 
further  simplify  both  the  analysis  and  the  computation  of  the  inversion,  we 
consider  a  constant  background  reference  medium.  By  examining  such  a 
rudimentary  problem,  we  are  able  to  focus  on  the  reliability  of  two- 
parameter,  constant  background.  Born  inversion  under  idealized  conditions. 
We  remark  also  that  this  thesis  is  organized  as  a  research  progression  which 
represents  the  evolution  of  our  understanding. 


DERIVATION  OF  THE  FORWARD  AND  INVERSE  EQUATIONS 


Introduction 

In  this  chapter,  we  derive  an  approziaate  technique  for  the  direct 
inversion  of  the  density  and  bulk  aodnlns  variations  of  a  stratified 
acoustic  earth.  As  mentioned  above,  we  reforaolate  the  method  of  Clayton 
and  Stolt  (1981)  to  accommodate  a  one-dimensional  medium.  We  therefore 
present  their  general  result  first.  In  their  paper,  Clayton  and  Stolt 
express  the  observed  wave  field  as  a  Lippmann-Schwinger  equation  which 
represents  the  sum  of  all  scattering  interactions  in  a  two-dimensional 
acoostic  medium.  The  resulting  integral  equation  is  equivalent  to  the 
Green's  theorem  expression  presented  by  Cohen  and  Bleiatein  (1979).  To 
invert  this  forward  equation,  Clayton  and  Stolt  assume  that  all  subsurface 
parameter  variations  are  small  and  apply  the  first  Born  approximation.  This 
procedure  establishes  an  approximate  linear  relationship  between  the  medium 
variations  and  the  observed  data,  and  it  is  identical  to  the  linearization 
inherent  in  the  method  of  Cohen  and  Bleistein.  Therefore,  although  the 
notation  of  this  derivation  is  similar  to  that  of  Clayton  and  Stolt,  the 
derivation  of  the  forward  equation  adheres  largely  to  the  work  of  Cohen  and 
Bleistein.  The  inversion  of  the  resulting  integral  equation  however  relies 
on  the  insights  of  Clayton  and  Stolt. 

Duriwntiom  of  the  Integral  Equation 

Wave  propagation  in  a  variable  density,  two-dimensional  acoostic  medium 
is  governed  by  the  equation: 


where  L  is  the  linear  wave  operator  defined  by  the  Biddle  expression  of  this 
dual  equality.  D  is  the  pressure  wave  field,  p  is  the  density,  K  is  the  bulk 
aodulus,  and  S(«)  is  the  Pourier  transfora  of  the  source-tiae  function.  In 
this  and  other  equations,  the  subscripts  g  and  s  denote  the  receiver  and  the 
source  locations,  respectively.  The  inversion  of  this  equation  for  the 
acoustic  aediom  paraaeters  is  in  general  a  nonlinear  problem.  To  overcome 
this  inherent  nonlinearity,  we  express  solntions  as  perturbations  about  a 
reference  medium  where  analytic  solutions  are  known.  The  equation 
representing  the  response  of  the  unperturbed  medium  is  given  by 


-6(x  -  x  )6(z  )S(id) 

g  s  s 


(2) 


where  is  the  linear  operator  in  the  reference  aediua,  Dr  is  the  field  in 
the  reference  aediua,  and  pf  and  are  the  reference  medium  parameters. 

Ve  now  rewrite  equation  (1)  so  as  to  introduce  the  scattering  potential 
V(x, z)  and  its  corresponding  perturbations  at  and  a*.  Becognizing  that  the 
scattering  potential  which  initiates  the  acattered  wave  field  is  aerely  the 
difference  of  the  linear  operators,  L  and  Lr,  we  add  and  subtract  Lr  to  L  in 
equation  (1).  After  collecting  teras,  we  find 
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(3) 
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» 


—  ▼  +  f-  +  V(x.z) 
pr  Kr 


D  =  -8(x  -  x  )6(z  )S (w) 

g  c  s 


where 


(4) 


and 

* 

<i)  a  a 

V(x.z)  -  -p-i  +  V  •  —  V  .  (5) 

K  p 

r  r 

Although  ax  and  aa  are  diaienalonl eas  parameters,  we  consider  them  as  the 

medium  variations  and  aeek  their  values  in  the  inversion  process.  To 
recover  the  actual  parameter  variations,  we  need  only  substitute  the 

inverted  values  of  ax  and  aa  into  equations  (4)  and  aolve  for  p  and  K. 

The  wave  field  repreaented  by  0  (recorded  at  the  geophone  location  xg 
due  to  the  aource  at  location  x#)  is  the  simple  integration  of  all 

scatterers  over  the  entire  subsurface  area.  This  total  wave  field  is  the 
sum  of  the  direct  wave  field  traveling  from  source  to  receiver  and  the  wave 
field  reflected  from  within  the  subsurface.  The  seismic  experiment  we  are 
considering  consists  only  of  the  measurement  of  the  scattered  or  reflected 
wave  field.  We  therefore  express  the  observed  field  or  the  data  as  the 

difference  of  the  total  field  and  the  field  in  the  reference  medium 

D(s  ,sg,»)  -  (D  -  IMSU)  (6) 

where  Ur  is  now  thought  of  aa  the  direct  wave.  The  preaense  of  S(a>) 
acknowledges  that  the  «'it  necessary  for  any  subsequent  inversion  are 
implicitly  band-limited. 

To  establiah  a  relationahip  between  the  data  and  the  parameter 


-  6  - 


variations,  we  subtract  equation  (2)  from  equation  (3) 

D  =  -V(Dr  +  D)S  (u>)  .  (7) 

In  order  to  provide  an  approximate  linear  relationship  between  the  surface 
observations  and  the  perturbations,  we  apply  the  first  Born  approximation 
and  replace  the  total  field  by  the  field  in  the  reference  medium.  As 
always,  the  Born  approximation  requires  that  the  perturbations  and  hence  the 
scattering  potential  remain  small.  Ve  are  justified  in  ignoring  the  product 
VO  appearing  on  the  right  hand  side  of  equation  (7)  as  it  is  of  lower  order 
than  the  product  VUr  (Bleistein,  1984).  Therefore,  the  linearized  equation 
is 


s 


V 


D  -  -VC  S (u) 
r 


(8) 


In  the  language  of  theoretical  physics,  the  above  procedure  is 
equivalent  to  truncating  the  Born  series  representation  of  equation  (7) 
after  the  second  term.  By  applying  this  linearization,  we  are  forcing  the 
inversion  to  model  all  observed  data  as  primary  reflections.  The  result  of 
this  is  that  the  inversion  cannot  properly  interpret  both  transmission 
losses  and  multiple  reflections.  The  addition  of  other  nonlinear  terms  of 
the  Born  series  would  include  such  phenomena  in  the  forward  equation,  but 


the  inherent  nonlinearity  woold  preclude  a  direct  inversion. 

It  is  convenient  to  express  equation  (8)  as  an  integral  equation  which 


incorporates  the  Green’s  function  representations  of  the  incident  and  the 
scattered  fields 


a  co 

D(xg,  x#,  <d)  ■J  “J  dz*  Gp(xg,0jx*,z'/w) 


•  V<x',z',w)G*(x’,z'|x#,0#w)S(w)  . 


(9) 


In  order  to  successfully  recover  the  components  of  reflectivity,  we 
have  considered  wave  motion  in  a  two-dimensional  medium.  In  keeping  with 
the  method  of  Clayton  and  Stolt  (1981),  we  consider  two-dimensional  Green’s 
functions  in  the  above  integral  equation.  The  first  Green's  function  in 
this  expression  represents  the  wave  propagating  from  the  scatterer  location 
to  the  receiver  location.  Similarly,  the  aecond  Green’s  operator 
corresponds  to  propagation  from  the  source  location  to  the  scatterer 
location.  Both  Green's  operators  represent  the  outgoing  solutions  in  the 
reference  medium  as  indicated  by  the  sign  and  are  derived  in  Appendix  A. 

According  to  equation  (9).  the  surface  wave  field  can  be  determined  if 
the  parameter  variations  aa  and  as  are  known.  Equation  (9)  is  therefore 
often  referred  to  as  a  forward  equation.  Our  goal  is  to  invert  this 
equation  so  that  the  parameter  varietions  are  expressed  in  terms  of  the 
known  data. 


Inversion  of  the  Integral  Bqnation  for  the  Redla  Paremeters 


For  our  specific  inversion  scheme,  we  ere  considering  only  constant 
background  reference  parameters,  pr  and  Ir»  This  approach  allows  a  direct 
inversion  of  equation  (9)  since  the  Green’s  functions  have  an  exact  analytic 


Substituting  equation 


D(x  > x  ,  to)  =  1/p  f  dx' 

S  *  r 


•  [  k’a^x'.z') 


(5)  into  equation  (9),  we  find 

dz  *  G  (x  ,  0  lx ',  z  '/«>) 
r  g  1 

+  V  •  aj(x,,z')V  1g*(x',  z  '  |x#*  0/w)S(w) 


(10) 


where  kf=  «/cr  .  Following  the  derivation  of  Clayton  and  Stolt.  we 
integrate  the  density  term  by  parts  to  obtain 
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(x*.  z* 

dz'  tk*ai(x', z')G^ 
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r  g 


(xg.0|x',z'/w)G*(x’,z'  |x#,0/u») 
•  VG*(x',  z '  |x#»  0/«)  )S  (w) 


(11) 


where  we  eliminate  the  boundary  terms  by  assuming  that  the  scattering 
potential  is  of  bounded  support.  Fourier  transforming  over  both  source  and 
receiver  locations  gives 


D(k  ,k  ,«)  * 
xg  xs 


1/p  J  dx'J  dz*  {k*ax(x’,  z ')Gr(kjg,  0  |x',  z */«)5r (x z '  |kxt,  Oaw) 


-  a.(x',z')Vff  (k  ,o|x'#z'/w)  •  ▼{?  (x’,z'|k  ,0/«))S(«)  (12) 

*  r  xg  1  r  1  is 


where  the  tilde  denotes  the  spatial  Fourier  transform. 

In  Appendix  A,  we  derive  the  analytic  expressions  for  the  Green's 


functions  represented  in  equation  (12). 


Inserting  these  results  into 


equation  (12)  and  applying  the  differential  operators  yields 


D(kxg'kx,'w) 


=  IT— C—  f  k*t  <k  -  k  ,-k  -  k  ) 

4kigkis  l*1**  xs  ig  zs 

+  (k  k  -  k  k  )t  (k  -  k  ,-k  -  k  )  ]s(w)  .  (13) 

zg  zs  xg  xs  *  xg  xs  zg  zs  J 


This  expression  establishes  an  approxiaate  linear  relationship  between  the 
Fourier  transform  of  the  scattered  wave  field  and  the  Fourier  transform  of 
the  aediua  variations.  As  pointed  out  by  Clayton  and  Stolt,  there  is  more 
than  enough  information  inherent  in  the  data  to  allow  an  inversion  for  the 
medium  parameters. 

In  this  investigation,  we  reduce  the  dimensionality  of  the  problem  by 
one.  That  is,  we  confine  our  attention  to  the  situation  where  the  medium 
parameters  depend  only  on  the  depth  variable.  As  a  result,  the  offset 
surface  reflection  data  must  be  collected  on  a  line-source,  line-receiver 
configuration,  and  only  one  shot  record  is  required.  Since  the  data  varies 
in  two  directions  while  the  medium  varies  in  only  one,  we  term  this  a  1.5 
dimensional  inversion.  In  Appendix  B,  equation  (13)  is  recast  in  terms  of 
the  one-dimensional  medium  parameters  to  give 


5(k  ,«)  «*  — l  \  k*I  (-2k  )  +  (2k*  -  k*)I  (-2k  )  Is (w)  .  (14) 

*  4k»  l  r  »  z  z  r  »  z  J 


Before  inverting  the  above  equation,  we  must  precondition  the  observed 
data.  If  the  data  were  a  recorded  time  section,  then  we  would  have  to 
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Fourier  transform  the  data  with  respect  to  x,  and  express  the  result  as  a 
function  of  k2 .  This  second  processing  step  would  be  achieved  by  applying 
the  dispersion  relation  and  resampling  the  resulting  function  for 
appropriate  values  of  •  Additionally#  we  would  also  have  to  deconvolve 
the  source  function  S(u>)  within  the  temporal  bandwidth  before  inverting. 
This  procedure  would  require  detailed  knowledge  of  the  source  function  which 
also  might  be  used  in  a  weighted  least  squares  solution.  Furthermore#  for  a 
fixed  angle  implementation  of  this  method  on  field  data,  the  deconvolution 
of  the  source  function  would  be  critical  for  accurate  inversion  results. 
For  the  purposes  of  this  investigation,  the  above  processing  steps  are 
unnecessary  since  we  are  considering  only  synthetic  data;  we  therefore 
express  our  synthetic  data  directly  as  a  function  of  k2  and  allow  all 
temporal  frequencies  to  be  of  equal  amplitude.  Thus,  the  forward  equation 
becomes 

D  (k  ,«o)  -  A  (k  ,w)»  (-2k  )  +  A  (k  ,w)» ,(-2k  )  (15) 

sz  »  z  1  z  »  z  *  z 


Since  the  medium  variations#  ax  and  sj»  independent  of  w,  we  may 

compute  Ax>  A,,  and  D$  for  two  distinct  values  of  «  and  solve  for  ax  and  a* 
directly.  With  temporal  frequency  as  the  free  parameter,  the  inversion 
relations  are  simplified  to  become: 
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a  a 
vs.  -  ii). 


(17a) 


(17b) 


By  expressing  the  inversion  relations  in  this  Banner,  ve  observe  the  near 
singnlar  behavior  of  these  functions  as  a>1  becomes  close  to  ua.  To  express 
Ij '  and  '  as  functions  of  depth,  we  inverse  Fourier  transform  equations 
(17). 

Since  the  choice  of  input  frequencies  is  somewhat  arbitrary,  we  could 
attempt  a  least-squares  or  an  average  least-squares  solution  over  frequency 
values.  This  approach  would  only  provide  improved  results  if  noisy 
synthetic  or  real  data  sets  were  being  inverted.  In  any  case,  the  least- 
squares  result  is  given  by  the  solution  to  the  equations: 


i }  v. 
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where  the  somations  are  taken  over  values  of  u.  The  derivation  of  the 


preceding  result  is  given  in  Appendix  C 


DERIVATION  OP  THE  SYNTHETIC  DATA 

Introduction 

In  this  chapter,  we  examine  the  synthetic  data  which  we  implement  to 
test  the  inversion,  le  choose  an  exact  analytic  form  for  the  data  so  as  to 
create  a  controlled  experiment.  This  control  not  only  provides  a  more 
confident  interpretation  of  subsequent  reconstructions,  but  also  it 
simplifies  the  mathematical  analysis.  To  conform  to  the  requirements  of  the 
inversion,  this  synthetic  data  is  computed  directly  in  the  Fourier  domain. 

Since  we  are  considering  a  one-dimensional  medinm,  the  synthetic  data 
is  constrained  by  a  line-source,  line-receiver  configuration.  Even  though 
the  medium  only  varies  with  depth,  the  data  varies  with  offset.  This 
feature  of  the  data  provides  the  information  necessary  for  the  separation  of 
the  components  of  impedance.  The  specific  form  of  this  data  is  derived  from 
the  point-source  results  given  by  Ewing,  Jardetzky,  and  Press  (1957).  In 
our  derivation,  we  follow  their  approach  and  develop  expressions  for  both 
single  and  double  interface  models.  In  summary,  the  resulting  analytic 
expressions  represent  the  two-dimensional  Fourier  transform  of  the  wave 
field  scattered  from  a  stratified  acoustic  medium  and  collected  at  many 
offsets . 


Derivation  of  the  Single  Interface  Date 

As  an  initial  attempt,  we  consider  an  analytic  expression  which 
represents  the  response  of  a  single  interface.  This  single  interface 


situation  is  modeled  by  joining  two  acoustic  constant  velocity  half-spaces. 


Since  we  are  inverting  for  one-dimensional  parameter  variations,  we  require 
only  one  shot  record/  the  experiment  providing  this  data  is  illustrated  in 
Fignre  1.  The  point  source  response  of  this  model  is  derived  by  Ewing, 
Jardetzky,  and  Press,  and  its  corresponding  line-source,  line-receiver 
expression  is  derived  in  Appendix  D: 
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w) 
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exp(-2Fih) 


(19) 


where  h  is  the  depth  of  the  reflector  below  the  source-receiver  line, 
6  =  P*/Pi»  end 


= 


The  right  side  of  equation  (19)  is  often  referred  to  as  the  plane-wave 

reflection  coefficient  for  the  single  interface  model  (Frisk,  1979),  and 
this  quantity  is  a  fnnction  of  both  the  incident  angle  and  the  input 

frequency. 

For  oar  pnrposes,  equation  (19)  must  be  written  as  a  function  of  kz. 
This  is  achieved  by  observing  the  dispersion  relation,  and  the  resulting 
expression  is 
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To  view  the  shot  record  associated  with  this  analytic  expression,  we 
have  only  to  apply  a  two-dimensional  FFT  to  equation  (19).  The  result  of 
this  procednre  for  a  particular  input  model  is  shown  in  Figure  2. 


Derivation  of  the  Doable  Interface  Data 


As  an  extension  of  this  single  interface  situation,  we  consider  a 
double  interface  model.  Although  this  model  includes  the  same  assumptions 
as  the  single  interface  case,  there  are  additional  features  which  may 
provide  insight.  Specifically,  the  two-interface  data  will  include  such 
nonlinear  effects  as  multiple  reflections  and  transmission  losses.  By 
inverting  this  forward  model,  we  should  be  able  to  examine  the  accuracy  with 
which  our  linear  inversion  reconstructs  nonlinear  data.  Furthermore,  as  the 
depth  and  complexity  of  the  probed  medium  increases.  Born  theory  predicts 
that  the  associated  error  will  increase.  Therefore,  by  examining  this  two- 
interface  model,  we  investigate  the  validity  of  our  linear  inversion  where 
it  is  predicted  to  breakdown. 

For  the  two-interface  forward  model,  we  consider  two  constant  velocity 
acoustic  half-spaces,  separated  by  a  constant  velocity  acoustic  layer.  As 
with  the  one-interfaci  model,  the  source  and  receivers  are  located  in  the 
upper  half-space.  The  geometry  for  this  model  is  illustrated  in  Figure  3. 
To  derive  the  analytic  expression,  we  follow  the  method  of  Ewing,  Jsrdetzky, 
and  Press  and  determine  the  displacement  potential  of  the  upgoing  waves  at 
the  receiver  depth.  This  is  accomplished  by  solving  the  system  of  equations 
which  satisfy  the  boundary  conditions  at  each  interface.  The  details  of 
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this  derivation  are  given  in  Appendix  E.  Thus,  the  synthetic  data  for  the 
two-interface  model  is  given  by  the  plane-wave  reflection  coefficient 
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81  +  R*  exp(-2^ha) 
1  +  RjBi  exp(-2^ha) 


exp(-2^lh1)  , 


(22) 


where  Rj  and  R2  »re  the  Rayleigh  reflection  coefficients  at  the  first  and 
second  interfaces,  respectively: 
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If  the  model  parameters  are  matched  op  across  either  interface, 
eqoation  (22)  reduces  to  equation  (19),  This  feature  allows  os  to  implement 
equation  (22)  for  either  the  single  or  the  double  interface  models. 

As  in  the  single  interface  model,  we  can  express  eqoation  (22)  as  a 


function  of  kzs  this  is  the  analytic  form  required  by  the  inversion 
relations  and  is  accomplished  by  rewriting  equations  (24)  as 
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A  = 


(26) 


To  a ccommoda te  the  Born  approximation,  we  may  choose  to  consider  a  form 
of  equation  (22)  which  excludes  multiple  reflections.  This  expression  is 
derived  in  Appendix  E  and  is  given  by 

B  (k  ,w)  =  R  exp(-2|f  h  )  +  R  Cl  -  R*)  expl-2(vh  +  tfh  >]  (27) 

where  (1  -  Rxs)  represents  the  transmission  loss  throngh  the  layer. 

An  example  of  a  shot  record  for  this  two-interface  model  is  shown  in 
Figure  4.  This  result  was  computed  by  applying  a  two-dimensional  FFT  to 
equation  (22). 

Before  inverting  equation  (15)  for  the  data  represented  by  equations 
(21)  and  (22),  we  note  that  the  synthetic  data  and  the  inversion  are  based 
on  two  different  Green's  function  conventions.  Equating  the  conventions 
used  by  Clayton  and  Stolt  (1981),  and  Ewing,  Jardetxky,  and  Press  (1957) 
gives  a  new  form  to  the  forward  equation: 
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A  (k  ,  to)  =  -2n  -  A  (k  ,  w) 
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(28) 


Thus,  to  invert  eqaation  (15)  for  the  synthetic  data  discussed  in  this 
chapter,  we  most  replace  eqaation  (16)  by  equation  (28), 
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AIT  IF  ACTS  IN  LINEAR  INVERSION 

Introduction 

Before  considering  examples  of  reconstr noted  acoustic  profiles,  we 
examine  the  source  of  artificial  parameter  variations  which  contribute  to 
two-parameter  reconstructions.  As  mentioned  above,  our  inversion  scheme  is 
based  on  a  linearization  of  the  nonlinear  relationship  between  the  model 
parameters  and  the  observed  data.  This  linearization  produces  errors  in  all 
subsequent  reconstructions.  For  instance,  a  variation  in  one  acoustic 
parameter  will  generate  a  false  variation  in  its  reconstructed  counterpart. 
We  term  these  artificial  parameter  changea,  inveraion  artifacts.  In  this 
chapter,  we  discuss  the  source  of  these  artifacts  in  the  context  of  linear 
inverse  theory.  Specifically,  we  outline  the  method  developed  by  Weglein, 
Violette.  and  Keho  (1985)  for  both  identifying  and  reducing  single 
interface  artifacts.  (In  fact,  this  thesis  represents  the  double  interface 
extention  of  their  work.)  In  this  chapter,  most  of  our  conclusions  are 
derived  from  the  analytic  form  of  the  synthetic  data.  We  begin  however  with 
a  general  analysis  of  the  differential  equations  relevant  to  our  method. 

Bxaaining  the  Differemtial  Equation 

Consider  the  differential  equation  governing  wave  motion  in  a 
stratified  acoustic  medium 
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Hie  spatially  transformed  wave  field  satisfies 
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and  can  be  rewritten  as 
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where  «o*/cr*  =  kx*  +  kz*  »  and 


V(z.ui) 


(32) 
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By  applying  the  Born  approximation,  we  assnae  that  the  total  field.  D,  can 
be  replaced  by  the  field  in  the  reference  medinm,  Br-  Since  the  field  in 
the  reference  medinm  satisfies 

0r(kx,  x,  w)  *  0  ,  (33) 

it  is  clear  that  equation  (31)  will  only  be  well  approximated  by  equation 
(33)  when 

k*ff  >>  VU,u)0  .  (34) 


For  the  fixed  frequency  implementation  outlined  in  the  previous  chapter,  kz 

f 

•  values  can  range  from  zero  to  «z/cr  where  tij  is  the  smaller  of  the  two 
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chosen  angular  frequencies.  Thus,  for  any  scattering  potential,  V,  there 
exists  kz  values  near  zero  for  which  equation  (34)  is  not  satisfied,  and  the 
Born  approximation  is  inappropriate.  Reconstructed  profiles  which  utilize 
low  kz  data  values  will  be  degraded  by  the  inherent  nonlinear  processes.  We 
should  therefore  eliminate  low  kz  data  values  prior  to  inversion.  In 
effect,  preconditioning  the  offset  dats  in  this  manner  allows  the  linear 
inversion  to  produce  more  exact  results. 


Analysis  of  Single  Interface  Inweraion  Artifacts 


In  the  previous  section,  we  examined  the  differential  equation  to 
determine  which  portion  of  the  data  in  the  wavenumber  domain  best  satisfies 
the  assumptions  in  our  linear  inversion.  We  found  that  even  though  no 
wavenumber  component  of  the  data  is  independent  of  nonlinear  effects, 
certain  portions  are  more  strongly  influenced  than  others.  To  better 
understand  this  observation,  we  now  consider  an  analytic  expression  for  the 
artifacts  generated  by  exact  synthetic  model  data.  As  a  simple  example,  we 
substitute  the  single  interface  result  given  by  equation  (21)  into  the 
inversion  relations  given  by  equations  (17). 

Consider  the  situation  where  only  the  bulk  modulus  varies  across  the 
interface.  In  this  case,  the  exact  value  of  the  density  variation  aa  is 
zero.  The  reconstructed  value  of  aa  is  determined  by  substituting 
into  equation  (17b)  with  px  *et  equal  to  pa: 
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exp(-2ikzh).  (35) 


This  relation  represents  a  density  change  at  z  =  h  and  is  identified  as  the 
source  of  the  density  artifact.  After  examining  this  equation,  we  conclude 
that  this  density  artifact  will  vanish  only  if  A  equals  zero.  Since  A  is 
zero  only  for  both  constant  density  and  bulk  modulus,  it  is  obvious  that  our 
inversion  scheme  will  always  predict  a  density  variation  for  the  case  of 
constant  density.  Alternatively.  it  can  be  shown  that  a  similar 
substitution  of  the  data  for  a  constant  balk  modulus  model  into  equation 
(17a)  leads  to  a  false  bulk  modulus  change.  Therefore,  even  though  the 
synthetic  data  has  an  exact  form,  the  reconstructed  acoustic  profiles 
contain  artificial  variations  which  are  traced  to  the  linearization  of  the 
inversion. 

To  understand  the  nature  of  these  artifacts,  we  examine  equation  (35) 
which  represents  the  artificial  density  variation  generated  by  a  bulk 
modulus  change.  Ve  notice  that  if  we  choose  tox  close  to  »*•  and  if 
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s  s 

kz  >>  to  A 


(36) 


then  the  density  artifact  diminishes/  that  is,  a,’  ->  0.  The  criterion 

expressed  by  eqaation  (36)  is  identical  to  equation  (34)  with  density  held 
constant.  As  before,  we  satisfy  this  requirement  by  eliminating  low  kz 
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offset  data  prior  to  inversion.  By  performing  a  similar  analysis  on  the 
bulk  modulus  artifact  generated  by  a  density  change,  we  confirm  that  the 
requirement  suggested  by  equation  (36)  will  indeed  reduce  this  artifact  as 
well.  We  conclude  that  the  condition  initiated  by  our  analysis  of  the 
differential  equation  is  verified  by  our  analysis  with  the  synthetic  data. 

If  we  eliminate  o>*  from  both  sides  of  equation  (36),  we  obtain 


cos*0  >> 


(37) 


where  0  represents  the  angle  of  incidence  with  respect  to  the  z-azis.  The 
condition  expressed  by  equation  (37)  suggests  that  the  data  corresponding  to 
large  0  are  responsible  for  the  inversion  artifacts.  Eliminating  low  k^ 
data  is  physically  equivalent  to  excluding  large  0  data  from  the  inversion. 
Note  also,  when  c*  <  c,,  equation  (37)  becomes  cos*0  >>  cos*0c  where  0C 
represents  the  critical  angle.  To  satisfy  this  condition,  we  must  insure 
that  0  <  @c  which  precludes  the  use  of  post-critical  data.  This  observation 
is  not  surprising  as  the  reflection  coefficient  becomes  imaginary  with 
modulus  equal  to  one  for  0  >  @c.  Equation  (37)  however  demands  more  than 
just  the  simple  elimination  of  inhomogeneous  waves  from  the  inversion.  To 
satisfy  equation  (37),  we  must  eliminate  the  incident  angles  which 
correspond  to  a  large  scattering  potential.  By  permitting  these  angles  to 
contribute  to  the  inversion,  we  are  allowing  the  scattering  potential  to 
become  large.  This  result  violates  the  primary  assumption  which  supports 
the  applicability  of  the  first  Born  approximstion.  Therefore,  to  diminish 
artificial  parameter  variations,  we  should  eliminate  the  data  which  are 
characterized  by  large  incident  angles  since  these  data  correspond  to  a 


t 
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significant  scattered  wave  field. 


In  his  two-parameter  Born  inversion  scheme,  Raz  (1981)  linearizes  the 
aconstic  reflection  coefficient  to  obtain  stable  results.  In  implementing 
Raz's  method,  Hanson  (1984)  remarks  that  the  reliability  of  this  method 
rests  largely  on  the  accuracy  of  this  approximation.  In  fact,  Hanson's 

results  are  favorable  only  when  the  variation  in  the  reflection  coefficient 
is  small.  At  large  angles,  the  reflection  coefficient  deviates  appreciably 
from  the  linearized  expression,  and  the  method  fails.  As  mentioned  by 
Clayton  and  Stolt,  the  amplitude  information  inherent  in  offset  data  allows 
the  proper  separation  of  the  components  of  impedance.  The  above  results 

however  suggest  that  the  large  6  information  can  seriously  degrade 

reconstructed  acoustic  profiles.  Ve  conclude  that  the  accuracy  of  these 
methods  relies  not  only  on  the  inherent  variation  in  the  reflection 

coefficient,  but  also  on  the  suitability  of  the  Born  approximation  in 
representing  this  offset  information.  Later,  we  will  consider  synthetic 
examples  which  illustrate  this  apparent  contradiction. 


Idemtif icatioa  of  Single  Interface  Inversion  Artifacts 


In  this  chapter,  we  have  presented  the  method  of  Veglein,  Violette  and 
Keho  (1985)  for  enhancing  the  inte<*pretabil  ity  of  reconstructed  acoustic 
profiles.  As  mentioned  above,  this  is  accomplished  by  preconditioning  the 
input  data  so  as  to  minimize  the  effects  of  nonlinear  processes.  In  their 
work,  Veglein,  Violette,  and  Eeho  also  discuss  a  method  for  identifying 
these  single  interface  artifacts/  we  outline  this  method  below.  Given  a 
data  set,  we  would  run  the  inversion  for  two  different  kz  bandwidths.  The 
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data  for  the  initial  inversion  would  include  all  permissible  k.  values  while 


the  second  inverted  data  set  would  be  filtered  to  eliminate  low  kz  values. 
Since  the  low  kz  data  is  highly  dependent  on  nonlinear  processes,  the 


broadband  reconstruction  would  be  less  reliable  than  the  band-limited 


FIXED  FREQUENCE  SINGLE  INTERFACE  RECONSTRUCTIONS 

Introduction 

In  this  chapter,  we  examine  the  reconstructed  acoustic  profiles  of  the 
single  interface  synthetic  data  discussed  in  Chapter  2.  Drawing  on  the 
analysis  of  the  previous  chapter,  we  implement  a  method  for  both  reducing 
and  identifying  artificial  parameter  variations.  In  all  these 
reconstructions,  the  reference  velocity  is  chosen  to  be  the  velocity  in  the 
first  medium.  Consequently,  the  Born-predicted  locations  of  the 
reconstructed  parameter  variations  are  always  correct,  and  the  Born-related 
artifacts  only  affect  the  amplitudes  of  the  reconstructions.  In  the 
subsequent  examples,  we  not  only  consider  the  deficiencies  of  our  linear 
theory,  but  also  we  emphasize  the  role  of  the  angularly  dependent  reflection 
coefficient  in  recovering  both  acoustic  parameters. 

Interpretation  of  the  Reeonatrnctiona 

Due  to  limitations  implicit  in  conventional  acquisition  techniques, 
seismic  data  is  band-limited  at  both  ends  of  the  spectrum.  This  loss  of 
information  impairs  all  subsurface  reconstructions  provided  by  direct 
inversion  procedures.  Specifically,  the  absense  of  high  spatial  frequencies 
limits  the  resolution  of  subsequent  images — structural  variation  within  the 
high  frequency  limit  is  not  detected.  An  equally  serious  limitation  is 
created  by  the  loss  of  low  frequency  information.  The  missing  low  spatial 
frequencies  contain  information  about  the  background  or  the  absolute  trend 
of  the  medium.  As  a  result  of  this  limitation,  all  reconstructed  parameter 
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changes  are  signaled  by  band-limited  step  functions  which  deviate  about  a 
constant  background  medium.  The  amplitudes  and  the  directions  of  actual 
parameter  changes  can  be  inferred  from  these  step  functions,  but  slowly 
varying  information  cannot  be  directly  recovered.  By  examining  the  impulse 


response 

of 

a  layered  earth 

model  for 

several  bandwidths. 

Hanson 

(1984) 

has 

demonstrated 

the  limitations 

implicit 

in  band-limited  data 

• 

For 

interpretat ional  purposes. 

Bleistein  and  Cohen 

(1982) 

choose 

to 

process 

the 

data  for  the 

derivative  of  the  variation 

ra  ther 

than 

the 

variation  itself.  This  result  arises  from  two  observations:  (1)  both  the 
amplitude  and  the  directon  of  a  parameter  change  are  more  easily  extracted 
from  a  band-limited  delta  function  than  from  a  band-limited  step  function, 
and  (2)  the  strength  of  a  parameter  change  can  be  directly  estimated  from 
its  corresponding  band-limited  delta  function  (Mager  and  Bleistein,  (1979)). 
In  this  investigation,  we  follow  the  suggestion  of  Bleistein  and  Cohen  and 
process  for  the  normal  derivative  of  the  perturbations.  Since  all 
reconstructions  are  contaminated  by  artifacts  of  the  linear  inversion,  we 
are  not  concerned  with  relating  these  reconstructed  amplitudes  to  the  actual 
input  parameter  changes.  We  note  however  that  the  peak  amplitudes  of  these 
sinc-like  functions  are  proportional  to  the  spatial  bandwidth 
(Nager  and  Bleistein,  (1979)).  In  presenting  the  reconstructions,  we 
normalize  the  reconstructed  profiles  to  the  largest  peak  value  of  ax  or  aa 
and  list  the  corresponding  maximum  values.  This  feature  aids  in  the 
comparison  of  different  reconstructions.  As  a  general  rule,  we  may  only 
compare  the  amplitudes  of  different  reconstructed  profiles  if  their  spatial 
bandwidths  are  equivalent. 
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Increasing  Velocity  Reconstructions 


In  this  section,  we  consider  the  inverted  ootpnt  of  three  different 
increasing  velocity  input  models.  These  models  simulate  :  (1)  constant 


density. 

(2) 

constant 

bulk  modulus. 

and 

(3)  constant 

normal 

inc idence 

impedance 

structures . 

Listed 

wi  th 

each 

reconstruction 

are 

the 

input 

velocities 

in 

feet  per 

se  cond, 

the 

input 

densities  in 

grams 

per 

cubic 

centimeter,  the  maximnm  incident  angle  in  degrees,  and  the  reflector  depth 
in  feet.  Also  included  is  a  plot  of  the  plane-wave  reflection  coefficient 
with  offset.  This  plot  not  only  demonstrates  the  nature  of  the  specific 
reflection  coefficient,  but  also  it  indicates  the  amount  of  offset  data 
included  in  the  inversion.  Note  that  when  kx  exceeds  the  smallest 
wavenumber  in  either  medium,  the  reflection  coefficient  becomes  imaginary. 
We  therefore  plot  the  reflection  coefficient  up  to  this  limit. 

In  these  examples  temporal  frequency  is  the  free  parameter,  and  the  two 
input  frequencies  are  listed  (in  Hz)  with  the  model  parameters.  The  smaller 
of  these  two  frequencies  determines  the  maximum  k  value:  k„_._  =  w,/c„. 
The  lower  limit  of  the  k2  bandwidth  is  fixed  by  the  maximum  offset  allowed 
in  the  inversion:  kzmin  =  krcos0Btx.  This  fact  follows  from  our  fixed 
frequency  formulation.  By  fixing  frequency,  we  must  sweep  through  a  range 
of  incident  angles  to  generate  a  suite  of  kz  values.  This  method  for 
establishing  a  spatial  bandwidth  has  its  shortcomings.  For  instance,  if  we 
reduce  offset  by  restricting  6,  then  we  necessarily  eliminate  some  low 
spatisl  frequencies  from  the  Inversion.  In  this  implementation,  the 
bandwidth  is  determined  by  both  the  choice  of  input  frequencies  and  the 
maximum  offset  included  in  the  inversion. 

To  conform  to  the  seismic  bandwidth,  we  generally  do  not  invert  for 
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frequencies  much  larger  than  60  Hz. 


Also,  we  choose  the  input  frequencies 


rather  close  together  since:  (1)  Born-related  artifacts  are  reduced,  and  (2) 
we  only  compute  the  inversion  for  kz  values  corresponding  to  the  smaller  of 
the  two  input  frequencies. 

As  a  first  model,  we  consider  a  10%  velocity  increase  with  density  held 
constant.  To  demonstrate  the  artificial  density  variation  produced  by  a 
bulk  modulus  change,  we  invert  this  model  with  all  pre-critical  offsets 
included  in  the  data.  This  reconstruction  is  shown  in  Figure  5.  As 
expected,  the  reconstructed  bulk  modulus  predicts  the  proper  parameter 
change  while  the  reconstructed  density  variation  records  a  substantial 
artifact.  This  artifact  is  a  plot  of  the  Fourier  transform  of  equation 
(35).  As  explained  in  the  previous  chapter,  the  origin  of  this  artifact  can 
be  traced  to  the  behavior  of  the  reflection  coefficient  at  far  offsets.  The 
plot  of  the  reflection  coefficient  for  the  data  inverted  in  Figure  5  details 
this  nonlinear  behavior.  To  reduce  this  artificial  density  variation,  we 
have  only  to  eliminate  the  far  offset  data  which  are  inconsistent  with  the 
Born  approximation.  The  reconstructed  profiles  for  this  case  are  presented 
in  Figure  6.  Notice  that  the  plot  of  the  reflection  coefficient 
demonstrates  the  reduction  of  offset  data  by  an  angle  of  10  degrees. 
Clearly,  the  density  artifact  has  diminished  relative  to  the  actual  bulk 
modulus  variation.  Ve  note  however  that  we  have  sacrificed  some  low  spatial 
frequencies  by  discarding  some  offset  information.  As  mentioned  above,  by 
reducing  the  spatial  bandwidth,  we  necessarily  reduce  the  amplitudes  of  the 
reconstructed  parameter  variations.  Therefore,  by  adhering  to  the 
condition  of  equation  (34),  ve  have  provider  an  inter pretabl e  acoustic 
profile  for  a  10%  velocity  change. 

The  preceeding  example  also  demonstrates  the  artifact  identification 
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procedure  of  Weglein,  Violette,  and  Keho  (1985).  By  comparing  the 
reconstructions  of  Figures  5  and  6,  we  observe  that  the  relative  amplitude 
of  the  density  variation  has  diminished  more  rapidly  than  the  bulk  modules 
variation.  Thos,  we  identify  the  density  change  as  an  artifact  of  the 
linear  inversion. 

As  a  second  example,  we  examine  the  inverted  profiles  of  an  input  model 
corresponding  to  a  10%  velocity  increase  for  constant  bulk  modulus:  c1*p1  = 
ci*P»-  To  yiew  the  associated  bulk  modulus  artifact,  we  run  the  inversion 
for  data  which  includes  all  pre-critical  offsets  (see  Figure  7).  As  before, 
we  reduce  the  Born-related  artifact  by  excluding  the  far  offset  (low  kz) 
data  from  the  inversion  (see  Figure  8).  Once  again,  we  acknowledge  that 
this  procedure  serves  to  diminish  the  relative  amplitude  of  the  artifact  and 
could  be  used  as  an  artifact  identification  procedure.  This  reconstruction 
however  has  a  different  character  from  the  previous  example.  By  comparing 
Figures  5  and  7,  we  notice  that  the  bulk,  modulus  artifact  does  not  dominate 
the  reconstruction  as  in  the  case  of  the  density  artifact.  This  result 
suggests  a  feature  of  two-parameter  Born  inversion  which  should  be  recalled 
when  interpreting  subsequent  reconstructions.  Simply  stated,  the 
reconstructed  density  profiles  are  less  reliable  than  the  reconstructed  bulk 
modulos  profiles.  Hanson  (1984)  has  also  documented  this  result  in  his 
implementation  of  Raz's  inverson  method.  Thus,  although  we  can  reduce 
artifacts  of  the  inversion,  we  recognize  that  both  parameter  reconstructions 
are  influenced  by  differing  levela  of  inaccuracy. 

As  pointed  out  by  Hinson,  the  most  accurately  inverted  profiles  are 
obtained  when  the  density  and  the  bulk  modulos  vary  in  a  manner  which 
reduces  impedance  contrasts.  This  situation  occurred  in  the  previous 
example.  As  an  additional  example,  we  consider  a  model  where  although  both 
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parameters  vary  across  the  interface!  the  normal  incidence  impedance  remains 


constant.  The  filtered  and  unfiltered  reconstructions  of  a  constant 
impedance  model  with  10%  parameter  changes  are  shown  in  Figures  9  and  10, 
respectively.  For  this  input  model,  the  expected  reconstruction  should 
consist  of  the  bulk  modulus  and  density  variations  moving  equal  amounts  in 
opposite  directions.  In  Figure  9,  we  notice  that,  once  again,  the  density 
variation  dominates  the  reconstruction.  As  before,  we  are  able  to  obtain  an 
interpretable  result  by  preconditioning  the  offset  data  prior  to  inversion 
(compare  Figures  9  and  10).  The  primary  reason  for  presenting  this  example 
is  to  demonstrate  the  usefulness  of  two-parameter  inversion  techniques.  By 
taking  advantage  of  impedance  variations  with  offset,  we  are  able  to 
separate  the  components  of  reflectivity.  Methods  which  utilize  coincident 
sour ce-re ce iv er  data  can  only  recover  normal  incidence  impedance  changes  and 
would  not  interpret  a  parameter  change  for  the  model  of  Figures  9  and  10. 


Deerensing  Velocity  Becoastractions 


When  the  velocity  decreases  across  the  interface,  there  are  no  post- 
critical  data  to  exclude  from  the  inversion.  This  allows  ns  to  test  the 
inversion  for  broadband  data.  As  an  initial  example,  we  consider  the  model 
of  Figure  5  with  the  velocities  interchanged.  This  situation  is  equivalent 
to  probing  the  model  of  Figure  5  from  below.  A  quick  look  at  the  plane-wave 
reflection  coefficient  for  this  model  reveals  that  its  magnitude  is 
identical  to  that  of  Figure  S.  In  Figure  11,  we  present  the  broadband 
reconstruction  of  this  model.  Clearly,  the  density  artifact  is  of  the  same 
order  as  the  actusl  bolk  modulus  variation.  Comparing  this  broadband 
artifact  with  its  counterpart  in  Figure  5,  we  observe  that  both  are  of  equal 
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size  relative  to  their  respective  bulk  modules  changes.  Since  the  spatial 
bandwidths  of  these  models  are  different,  this  equivalence  might  be  traced 
to  the  equivalence  of  the  reflection  coefficients.  To  test  this,  we  invert 

the  model  of  Figure  12  for  data  which  correspond  to  the  reflection 

coefficient  of  Figure  6.  In  comparing  the  relative  sizes  of  the  resulting 
density  artifacts,  we  conclude  that  they  indeed  are  of  equal  amplitude. 
Therefore,  the  amplitudes  of  inversion  artifacts  are  determined  by  the 

nonlinear  nature  of  the  reflection  coefficient.  By  adhering  to  the 

criterion  suggested  by  equation  (37),  we  can  reduce  these  deleterious 
effects. 


Constant  Velocity  Ceconstrnctions 

As  a  final  example  which  illustrates  the  usefulness  of  two-parameter 
inversion  techniques,  we  consider  a  constant  velocity  model.  For  this 
situation,  the  reflection  coefficient  is  a  constant  which  does  not  vary  with 
offset.  Ve  therefore  gain  nothing  by  filtering  the  offset  data  prior  to 
inversion.  This  example  provides  the  only  situation  where  our  approximate 
inversion  agrees  with  the  exact  result.  That  is,  since  there  is  no  velocity 
constrast,  the  Born  approximation  is  exactly  satisfied  (\x  =  0)  .  As  a 
result,  both  tne  Born  and  exact  theories  predict  that  both  parameter 
variations  have  equal  amplitude  and  polarity.  An  example  of  a  broadband 
constant  velocity  reconstruction  is  illustrated  in  Figure  13.  Clearly,  this 
two-parameter  inversion  scheme  can  accurately  predict  these  subsurface 
parameter  changes. 
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Least-Squares  Reconstructions 


As  mentioned  above,  inversion  results  improve  if  the  input  frequencies 
are  chosen  close  together.  There  is  however  no  criterion  which  governs  the 
specific  choice  of  these  frequencies.  In  the  previous  examples,  we  chose  60 
and  62  Hz,  because  these  values  both  conform  to  the  seismic  bandwidth,  and 
provide  reasonable  inversion  results.  To  demonstrate  that  the  specific 
choice  of  frequencies  is  arbitrary,  we  consider  a  leaat-sqoares  solution 
over  a  range  of  frequency  values.  In  Figure  14,  we  present  a  least-squares 
reconstruction  of  the  model  shown  in  Figure  6.  These  results  are  obtained 
by  inverting  equation  (18)  with  the  summation  extending  from  55  to  65  Hz. 
In  comparing  Figures  6  and  14,  we  acknowledge  that  the  reconstructions  are 
similar  in  character  and  would  lead  to  equivalent  interpretations.  Their 
only  difference  lies  in  the  frequency  content  of  the  reconstructed  profiles. 
The  similarity  of  these  results  is  not  alarming.  Ve  do  not  expect  the 
least-squares  procedure  to  improve  the  inversion  of  synthetic  data  since 
this  method  is  only  useful  when  considering  noisy  synthetic  or  field  data 
sets.  The  above  example  merely  demonstrates  that  inversion  results  are  not 
strictly  dependent  on  the  choice  of  input  frequencies.  As  mentioned  above, 
the  two  frequencies  should  be  chosen  close  together,  and  should  be  large 
enough  to  provide  a  sufficient  bandwidth  for  stable  computation. 


FIXED  FIBQOENCI  DOUBLE  INTERFACE  RECONSTRUCTIONS 


Introduction 

In  all  wave  equation  methods  which  transform  sorface  reflection  data 
into  a  map  of  subsurface  reflectors,  the  selection  of  a  background  velocity 
profile  determines  the  accuracy  of  the  results.  For  Born  inverse  methods, 
this  reference  velocity  affects  not  only  the  migration  of  the  image,  but 
also  the  recovery  of  the  amplitude  information.  These  results  are 
demonstrated  in  a  simple  example  by  Weglein  and  Gray  (1983).  In  this 
example,  they  consider  the  single  parameter  Born  reconstruction  of  a  one¬ 
dimensional  medium  for  various  constant  background  velocities.  They 
conclude  that  constant  background  Born  techniques  cannot  simultaneously 
determine  both  the  location  and  the  size  of  a  single  parameter  change. 

As  mentioned  above,  our  inversion  scheme  assumes  that  the  velocity  in 
the  first  medium  is  the  constant  background  velocity.  Consequently,  the 
locations  of  both  parameter  variations  are  correct  for  the  first  reflector 
while  the  amplitudes  are  necessarily  incorrect.  In  addition,  these 
amplitudes  are  further  complicated  by  the  apparent  link  between  the 
parameter  variations.  The  deficiencies  of  Born  methods  do  not  improve  with 
depth.  In  fact,  both  the  locations  and  the  amplitudes  of  the  parameter 
variations  for  the  second  reflector  will  be  incorrect.  In  this  chapter,  we 
examine  the  reconstructions  of  the  double  interface  synthetic  data  derived 


in  Chapter  2. 

When  examining 

these  reconstructions,  we 

use 

the  first 

reconstructed 

reflector  a*  a 

measure  for  the  accuracy 

of 

the  second 

reconstructed 

reflector.  For 

interpretations!  purposes. 

it 

is  safe  to 

assume  that  the  amplitude  of  the  first  change  is  correct  if  the  data  have 


been  filtered.  The  following  examples  will  expose  the  deficiencies  inherent 
in  the  formulation  of  onr  approximate  inversion.  As  before,  we  will  attempt 
to  rectify  inversion  artifacts  by  filtering  the  inpot  data  prior  to 
inversion. 


Examples  of  Doable  Interface  ieconstroctiona 

To  conform  to  realistic  earth  models,  we  considered  10%  parameter 
changes  for  the  single  interface  model.  By  filtering  the  input  data,  were 
able  to  obtain  interpretable  reconstroctions .  For  the  double  interface 
model,  we  begin  by  considering  two  consecutive  10%  parameter  changes.  In 
Figure  15,  we  present  the  reconstruction  of  a  constant  density  model  with 
the  layer  separating  two  half-spaces  of  equal  but  lower  velocity.  Clearly, 
the  inversion  is  unstable  for  this  input  model/  the  density  artifact 
associated  with  the  second  interface  dominates  the  reconstruction.  To 
improve  this  reconstruction,  we  adhere  to  the  previous  strategy  and 
eliminate  far  offset  data  prior  to  inversion  (tee  Figure  16).  In  comparing 
Figures  15  and  16,  we  observe  that  this  procedure  does  not  improve  the 
interpretability  of  the  reconstruction— eliminating  low  kz  data  merely 
reduces  the  spatial  bandwidth.  Furthermore,  the  artifact  identification 
procedure  outlined  in  Chapter  3  ia  not  effective  in  delineating  the  double 
interface  density  artifact.  We  conclude  that  this  procedure  is  limited  to 
situations  where  the  inverted  wave  field  is  accurately  modeled  above  the 
reflector  of  interest.  Ve  note  also  that  the  density  artifact  incorrectly 
establishes  the  location  of  the  second  interface.  Born  theory  however  would 
predict  a  location  shallower  than  3000  feet  rather  than  deeper.  Ihns,  it 
appears  that  this  inversion  scheme  produces  unstable  results  for  this  model. 
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and  these  results  are  not  improved  by  further  filtering  the  inpot  data. 

Since  this  method  is  based  on  a  weak  scatterer  approximation,  we  expect 
stable  reconstr actions  for  smaller  parameter  changes.  In  Fignres  17  and  18, 
we  present  the  reconstructed  profiles  of  the  model  of  Figure  IS  with  the  10% 
velocity  changes  replaced  by  1%  and  0.1%  changes,  respectively.  In  all 
three  models,  the  maximum  offset,  the  input  frequencies,  and  the  reflector 
depths  are  equivalent.  Tnrning  first  to  Figure  17,  we  observe  that  although 
the  density  artifact  continues  to  dominate  the  reconstruction,  the  overall 
reaults  are  more  atable.  Note  that  the  reflection  coefficient  is  trivial 
for  all  offsets  in  this  inversion,  yet  the  reconstruction  is  inaccurate. 
This  result  further  demonstrates  that  the  low  filter  has  a  limited 

application  for  the  two-interface  situation.  When  the  parameter  changes  are 
reduced  by  another  order  of  magnitude  (see  Figure  18),  the  reconstruction  is 
interpretabl e/  however,  the  density  artifact  located  at  the  second  interface 
still  exists.  For  this  model,  the  reflection  coefficient  is  too  small  to  be 
distinguished  from  the  kx~axis.  From  the  proceeding  examples,  we  conclude 
that  this  inversion  method  cannot  properly  interpret  two  successive  real 
earth  parameter  changes.  As  shown,  the  inversion  results  are  only  accurate 
for  a  cumulative  change  of  0.2%. 

In  the  previous  chapter,  we  examined  a  constant  bulk  modulus  model  to 
illustrate  that  inversion  results  improve  for  parameter  changes  which 
minimize  impedance  contrasts.  In  Figures  19  through  21,  we  present  the 
reconstructed  profiles  of  double  interface,  constant  bulk  modulus,  models. 
These  models  are  consecutively  equivalent  to  the  models  of  Figures  15,  17, 
and  18  with  only  the  densities  adjusted  to  simulate  constant  bulk  modulus. 
For  example,  the  spatial  bandwidth,  the  velocity  contrasts,  and  the 
reflector  locations  are  equal  for  the  models  of  Figures  15  and  19/  however. 
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the  inpnt  density  profile  of  Figure  19  produces  ■  constant  bulk  modulus.  By 


equating  all  but  the  input  densities,  we  insure  that  the  limitations  imposed 
by  the  Born  approximation  are  identical.  This  allows  a  aimple  comparison 
between  the  two  models.  Similarly,  the  constant  denaity  models  of  Figures 


17  and  18  can  be 

compared  with 

the  conatant  bulk  modulus  models 

of 

B 

Figures  20  and  21, 

reapect ively . 

Thus,  in  considering 

these 

six 

reconstructions,  we 

can  compare  results  within  each  constant 

parameter 

set 

or  against  its  counterpart  in  the  other  set. 

_  Considering  Figure  19  first,  we  notice  that  the  density  change  located 

at  the  second  interface  dominates  the  reconstruction.  This  change  should  be 
of  equal  amplitude  and  opposite  polarity  as  the  change  which  delineates  the 

I 

first  interface.  We  therefore  consider  this  reconstruction  to  be  unstable. 
In  fact  the  character  of  this  reconstruction  resembles  the  unstable  constant 
density  reconstruction  of  Figure  IS.  The  differing  density  profile  only 
fi  affects  the  response  at  the  first  interface.  We  conclude  that  for  velocity 

contrasts  of  this  magnitude,  the  inversion  results  are  unstable,  regardless 
of  the  input  density  profile.  It  is  also  important  to  note  that  the 

flj  reconstructed  density  is  responsible  for  the  instability.  This  result 

reinforces  the  previous  statement  that  the  reconstructed  density  is  less 
reliable  than  the  reconstructed  bulk  modulus. 
w  Next,  we  examine  the  conatant  bulk  modulua  reconstruction  of  Figure  20. 

Although  the  velocity  structure  in  this  model  corresponds  to  the  constant 
density  model  of  Figure  17,  this  reconatruction  is  superior.  In  Figure  17, 
the  density  artifact  at  the  second  interface  dominates  the  reconstruction. 
In  Figure  20,  the  density  change  at  the  second  interface  is  only  slightly 

incorrect  in  amplitude  and  location.  Alao,  the  bulk  modulus  artifact  is 

A  quite  small.  Thus,  for  1%  velocity  changes,  the  input  density  profile  of 


\ 
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Figure  20  reduces  the  impedance  contrasts  and  improves  the  inversion.  In 
fact,  the  accuracy  of  this  reconstruction  is  comparable  to  that  of  Figure 
18,  even  though  the  velocity  contrasts  are  ten  times  greater.  This  further 
suggests  that  artificial  bulk  modulus  variations  are  less  damaging  than 
artificial  density  variations. 

Finially,  we  consider  the  reconstrnction  of  Figure  21.  As  expected, 
this  result  is  almost  exact/  the  bnlk  modulus  artifact  is  trivial.  In 
comparing  all  three  constant  bulk  modulus  reconstroctions  with  all  three 
constant  density  reconstroctions,  we  conclude  that  the  former  set  are  more 
accurate  than  the  latter.  This  result  is  attributed  both  to  the  reduction 
of  impedance  across  the  interfaces,  and  the  controlling  influence  of  the 
density  variation  in  all  reconstroctions.  Ve  most  admit  however  that  the 
inversion  is  unstable  for  the  realistic  earth  models  of  Figures  15  and  19. 

Returning  to  the  constant  density  model  of  Figure  15,  we  consider  two 
equivalent  models  but  for  the  depth  of  the  second  interface.  In  Figure  22, 


we  present 

the 

reconstruction  of 

this  model  with 

a 

shallower 

second 

refl ector . 

The 

model  of 

Figure  23 

is  characterized 

by 

a  deeper 

second 

reflector. 

In 

comparing 

these  three 

reconstructions. 

we 

conclode 

the  the 

error  associated  with  the  approximate  reconstrnction  of  the  first  interface 
increases  with  depth.  Therefore,  not  only  is  this  fixed  frequency  procedure 
unstable  for  two  realistic  parameter  changes,  but  also  the  accuracy  of  all 
reconstructions  depends  on  the  depth  of  the  second  reflector.  Instead  of 
speculating  as  to  the  origin  of  these  deficiencies,  we  take  advantage  of  the 
exact  form  of  the  model  data  to  analyze  inversion  artifacts.  Before  ve  do 
this,  we  reformulate  the  inversion  to  both  improve  stability  and  simplify 
analysis . 
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FIXED  ANCLE  INVERSION 


Introdact ion 


In  the  previous  inversion  scheme,  temponl  frequency  is  chosen  ss  the 
free  variable  when  solving  for  the  parameter  variations.  In  that 
formulation,  the  necessary  suite  of  values  is  generated  by  sweeping 

through  a  range  of  incident  angles  for  each  fixed  wave  vector.  A 
prohibitive  consequence  of  this  procedure  is  that  the  elimination  of  the 
detrimental  far  offsets  severely  reduces  the  spatial  bandwidth.  As  a 
result,  the  inversion  becomes  unstable.  To  overcome  this  instability,  we 
reformulate  the  inversion  in  terms  of  fixed  incident  angles.  Instead  of 
choosing  two  fixed  frequencies,  we  now  choose  two  fixed  incident  angles  when 
solving  for  the  parameter  variations.  The  kz  values  are  computed  by 
establishing  a  temporal  frequency  bandwidth  for  the  fixed  angle  wave  vector. 
As  we  will  see.  this  formulation  achieves  greater  stability  as  the  offset 
information  and  the  spatial  bandwidth  are  no  longer  linked.  In  this  new 
formulation,  we  can  choose  incident  angles  which  correspond  to  Born- 
compatible  reflection  coefficients  with  no  restriction  on  spatial  bandwidth. 
This  procedure  is  not  unique.  In  his  time  domain  formulation,  Raz  (1981) 
utilizes  two  distinct  offsets  when  recovering  both  velocity  and  density. 
Therefore,  in  this  chapter,  we  reformulate  the  inversion  to  improve  the 
stability  of  the  output.  To  demonstrate  this  improvement,  we  re-invert  the 
models,  previously  considered  in  the  fixed  frequency  formulation.  As  an 
additional  feature,  this  reformulation  simplifies  the  computation  of  the 


inverted  profiles, 
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Reformulation  of  the  Inversion 


In  this  section,  we  rewrite  the  forward  equation  in  terms  of  the 
incident  angle  9.  This  is  achieved  by  substituting  the  definition 

k  =  — —  cos0  (3  8) 

z  c 

r 


in  equations  ( 28)  : 


A^e) 


-n 

cos*0 


Aa ( 0)  =  -2n  -  Ai(0) 


(39) 


where 


Ds(kz.0) 


A  (0)f  (-2k  )  +  A  (0)«r  (-2k  ) 
i  i  z  a  *  z 


(40) 


To  invert  equation  (40)  for  the  transformed  medium  variations,  we  compute 
At,  A,,  and  Dg  for  two  values  of  0  and  solve  for  tx  and 


a ' ( -2k  )  =  t'(-2k  )  +  - 
»  z  *  z  n 


D  (k  ,0  )  -  D  (k  ,0  ) 
a  z  a  a  z  a 


sec 


ae c  9, 


(41a) 


t'(-2k  ) 
*  z 


1_ 

2n 


sec 


*0  D  (k  ,0  )  -  aec*0  5  (k  ,0  ) 


a  a  z  a 


a  a  z  a 


a  a 

sec  9X  -  sec  0, 


(41b) 
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Clearly,  this  inversion  procedure  is  computationally  simpler  than  the 
previous  method  since  equations  (39)  are  independent  of  k2 .  For  a  specific 
6,  Aj  and  Aa  are  numbers  which  weight  each  spatial  component  of  the  data 
equally.  This  result  not  only  simplifies  the  computation,  but  also  it 
reduces  the  complexity  of  future  mathematical  analysis.  Once  again,  we  note 
the  near  singular  behavior  of  the  denominator  of  equations  (41)  for  61  close 
to  6a.  As  before,  we  have  only  to  inverse  Fourier  transform  equations  (41) 
to  express  the  parameter  variations  as  a  function  of  depth. 

To  compute  the  synthetic  data  for  this  method,  we  generate  the  kz  I 

values  by  fixing  the  incident  angle  and  sweeping  through  a  specified  range 
of  temporal  frequencies.  For  this  fixed  angle  formulation,  we  re-express 
the  single  interface  data  of  equation  (21)  as 


Ds(kz,G) 


t 

c  A 
r  x 

COS  *0 


cX 

cos  0 


exp(-2ik  h)  , 
r 


and  the  double  interface  data  of  equation  (22)  as 


(42) 


6  (k  .0) 
s  x 


Bx(0)  +  8j(0)exp(-2*fjhi) 
1  +  B1(G)8J(6)exp(-2*»lh1) 


exp(-2|/Jh1) 


(43) 


where  Bt  and  B,  are  the  Bayleigh  reflection  coefficients  at  the  first  and 
aecond  interfaces  respectively,  and 


cos  6 


From  these  expressions,  we  observe  that  the  Rayleigh  reflection  coefficients 


■re  also  independent  of  k2, 


Examples  of  Fixed  Angle  Reconstructions 


In  this  section,  we  examine  fixed  angle  reconstructions  of  single  and 
double  interface  models  presented  in  Chapters  4  and  5.  To  compare  the  two 
inversion  procedures,  we  equate  both  the  input  model  parameters  and  the 
spatial  bandwidths  before  inverting.  In  all  of  these  examples,  we  examine 
the  reconstruction  of  a  constant  density  profile/  this  allows  ns  to  test  the 
inversion  for  unfavorable  conditions.  As  an  initial  example,  we  consider 
the  fixed  angle  reconstruction  of  the  single  interface  model  shown  in  Figure 
6  (see  Figure  24).  Note  that  we  now  list  the  input  angles  (in  degrees)  and 
the  temporal  bandwidth,  instead  of  the  input  frequencies  and  the  maximum 
offset.  In  comparing  Figures  6  and  24,  we  observe  that  the  bulk  modulus 
variations  are  strikingly  similar  while  the  density  artifact  is  smaller  in 
the  fixed  angle  reconstruction.  Ihis  improvement  results  from  our  ability 
to  confine  the  wave  vector  to  a  region  which  corresponds  to  a  small 
reflected  wave  field.  In  other  words,  in  the  fixed  angle  implementation,  we 
are  able  to  satisfy  equation  (37)  by  choosing  incident  angles  close  to  zero. 
Consequently,  the  reflection  coefficient  satisfies  the  Born  approximation. 
This  idea  can  be  seen  by  comparing  the  plots  of  the  reflection  coefficients 


in  Figures  6  and  24.  In  the  fixed  angle  implementation,  the  reflection 


coefficient  maintains  the  low  values  determined  by  the  incident  angle  of  10 
degrees,  whereas  its  counterpart  in  Figure  6  traces  oot  values  corresponding 
to  55  degrees.  Clearly,  the  plot  of  the  reflection  coefficient  for  the 
fixed  angle  implementation  appears  linear  (constant).  The  selection  of 
incident  angles  in  the  fixed  angle  method  is  arbitrary.  Ve  note  however 
that  as  in  the  fixed  frequency  formulation,  inversion  artifacts  diminish  if 
the  input  parameters  (the  incident  angles  in  this  case)  are  chosen  close  to 
each  other.  In  all  subsequent  reconstructions,  we  select  the  angles  of  5 
and  10  degrees  to  comply  with  the  Born  approximation.  Therefore,  the 
flexibility  of  the  fixed  angle  procedure  provides  more  accurate  results  as 
the  spatial  bandwidth  is  no  longer  linked  to  the  offset  information. 

As  demonstrated  in  the  previous  example,  the  fixed  angle  implementation 
improves  single  interface  reconstructions  by  farther  discriminating  against 
nonlinear  portions  of  the  offset  data.  We  now  apply  this  method  to  invert 
doable  interface  models  and  compare  these  resnl ts  with  their  corresponding 
fixed  frequency  reconstructions.  As  a  first  example,  we  invert  the  constant 
density  model  of  Figure  15  by  the  fixed  angle  method/  this  result  is  shown 
in  Figure  25.  As  before,  we  are  able  to  compare  these  reconstructions  since 
both  inversion  methods  are  computed  for  equivalent  spatial  bandwidtbs.  In 
comparing  Figures  15  and  25,  we  notice  that  the  density  artifact  at  the 
second  interface  dominates  both  reconstructions.  However,  the  instability 
inherent  in  the  fixed  frequency  reconstruction  is  not  present  in  the  fixed 
angle  reconstruction.  In  the  fixed  angle  reconstruction,  the  Born-predicted 
location  of  the  density  artifact  is  correct.  Thus,  although  the  fixed  angle 
technique  enhances  the  character  of  this  reconstruction,  the  results 
indicate  that  Born  methods  cannot  successfully  reconstruct  both  acoustic 
parameters  for  10%  velocity  changes. 
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To  continue  this  comparison,  we  consider  two  other  constant  density 
models  with  1%  and  0.1%  velocity  changes,  respectively.  In  Figure  26,  we 
present  the  fixed  angle  reconstruction  of  the  1%  velocity  contrast  model  of 
Figure  17.  In  comparing  these  two  reconstructed  profiles,  we  observe  that 
the  fixed  angle  profile  is  more  accurate  than  the  fixed  frequency  profile. 
In  the  fixed  angle  reconstruction,  the  density  artifact  located  at  the 
second  reflector  is  one  third  the  size  of  its  fixed  frequency  counterpart. 
Nevertheless,  this  artifact  is  still  appreciable.  To  complete  this  series 
of  comparisons,  we  consider  the  fixed  angle  reconstruction  of  the  0.1% 
velocity  contrast  model  shown  in  Figure  18  (see  Figure  27).  Although  both 
reconstructed  profiles  are  interpretable  in  terms  of  the  input  model,  the 
results  of  Figure  27  are  superior.  As  before,  the  density  artifact  is 
smaller  for  the  fixed  angle  reconstruction.  Therefore,  even  though  the 
aforementioned  advantages  of  the  fixed  angle  method  improve  the 
reconstructions,  the  inversion  method  is  still  unreliable  for  real  earth 
parameter  changes. 

In  the  fixed  frequency  doable  interface  examples  of  the  previous 
chapter,  we  presented  evidence  which  suggested  that  errors  associated  with 
the  incomplete  (Born)  reconstruction  of  the  first  interface  propagate  with 
depth.  To  verify  this  statement,  we  consider  fixed  angle  reconstructions  of 
a  constant  density  input  model  for  various  depths  of  the  second  reflector. 
We  choose  the  input  model  of  Figure  26.  Note  that  in  the  fixed  frequency 
demonstration  of  this  phenomenon,  we  considered  0.1%  velocity  contrasts. 
For  the  fixed  angle  example,  we  consider  1%  velocity  contrasts  across  the 
interfaces.  In  Figures  28  and  29,  we  present  the  reconstructed  profiles  of 
the  model  of  Figure  26  with  the  second  reflector  located  at  2000  and  4000 


feet,  respectively.  In  comparing  these  three  figures,  we  conclude  that  the 
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error  does  in  fact  increase  with  depth.  This  feature  of  Born  inversion  is 
equally  damaging  for  both  inversion  procedcres. 

In  comparing  the  fixed  angle  reconstructions  of  Figores  26.  28,  and  29, 
we  observe  that  the  error  affects  not  only  the  amplitude  but  also  the 
location  of  the  parameter  changes  associated  with  the  second  interface. 
This  result  follows  from  the  assumption  of  a  constant  background  velocity. 
By  assuming  a  constant  background  velocity,  we  force  the  inversion  to  model 
all  waves  as  if  they  are  propagating  at  the  reference  velocity.  In  this 
example,  the  reference  velocity  is  the  velocity  in  the  first  medium.  Since 
this  velocity  is  slower  than  the  velocity  in  the  second  medium,  the 
reconstructed  parameter  changes  are  located  above  their  expected  locations. 
This  featore  of  the  inversion  is  more  obvious  in  the  fixed  angle  example 
than  in  the  corresponding  fixed  frequency  example  because  the  velocity 
contrasts  are  larger.  In  the  next  chapter,  we  consider  the  origin  of 
inversion  artifacts  in  the  context  of  the  preceeding  discussion. 

In  the  above  examples,  we  equate  spatial  bandwidths  between  inversion 
procedures  to  allow  a  direct  comparison  of  the  amplitudes  of  the 
reconstructed  profiles.  This  apparent  limitation  is  in  opposition  to  the 
flexibility  of  the  fixed  angle  approach.  In  the  fixed  angle  approach,  we 
may  choose  any  spatial  bandwidth  for  any  incident  angle.  As  a  result,  we 
are  free  to  aelect  input  data  in  accordance  with  the  assumptions  implicit  in 
the  linear  inversion.  Unlike  the  fixed  frequency  method,  the  stability  of 
the  results  are  not  compromised  by  the  bandwidth  restrictions  effected  by 
the  Born  approximation. 
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INVERSION  ARTIFACTS  REVISITED 


Introduction 

In  Chapter  3,  we  examined  the  analytic  form  of  the  artificial  parameter 
variations  associated  with  the  Born  reconstr action  of  the  first  interface. 
Ve  showed  that  these  artifacts  conld  be  rednced  by  emphasizing  certain 
portions  of  the  offset  data  in  the  inversion  process.  In  fact,  these  single 
interface  artifacts  can  practically  be  eliminated  through  the  flexibility  of 
the  fixed  angle  implementation.  In  addition  to  these  false  parameter 
variations.  Born  theory  inaccnrately  predicts  the  size  of  the  actnal 
parameter  change  at  the  first  interface,  even  if  the  velocity  in  the  first 
mediom  is  selected  as  the  reference  velocity.  Fortunately,  this  amplitude 
error  is  trivial.  When  considering  the  double  interface  reconstructions,  we 
demonstrated  that  the  errors  associated  with  the  incomplete  reconstruction 
of  the  first  interface  increase  with  depth.  Thus,  even  though  these  errors 
are  not  severe  at  the  first  interface,  they  are  significant  at  the  second 
interface.  In  this  chapter,  we  examine  the  analytic  expressions  of  the 
artifacts  which  characterize  the  reconstructed  changes  at  the  second 
interface.  This  analysis  both  verifies  previous  conclusions  and  provides 
insight  into  a  possible  correction  of  errors  implicit  in  two-parameter, 
constant  background.  Born  inversion. 

Analysis  of  Doable  Interface  Artifacts 

As  in  the  single  interface  analysis,  we  substitute  constant  parameter 
model  data  into  the  corresponding  inversion  relation  to  obtain  an  expression 


for  the  artifact.  Specifically,  we  will  again  consider  the  density  artifact 
prodoced  by  a  balk  model  os  change.  For  this  doable  interface  analysis,  we 
exploit  the  simplicity  offered  by  the  fixed  angle  implementation.  Since  the 
multiple  reflections  do  not  interfere  with  the  primary  reflections,  we 
otilize  the  fixed  angle  expression  of  the  model  data  given  by  equation  (27) 

B^k^.0)  =  Bj(0)exp(-2yihl)  +  BJ(0)T(0)exp[-2(|(1h1  +  i^h,)]  (45) 

where  ffx  md  j/a  are  computed  according  to  equations  (44)  and  I  represents 
the  transmission  loss  through  the  layer.  In  this  expression,  each  term 
represents  the  primary  reflection  event  for  each  interface.  Substituting 
equation  (45)  with  a  constant  density  profile  and  an  arbitrary  reference 
velocity  into  equation  (41b),  we  obtain  an  expression  for  the  density 
artifact : 


'  1 
*  <~2k  )  =  - — - - - - 

z  2n(sec  ©i  -  sec  0#) 

•  sec  exp[-2iki(l  +  c^A/cos  &x)  *  h4J 

-  sec  0  B  (0  )  exp[-2ik  (1  +  c  A/cos*0  )  ^h  ] 

ii>  2  r  *  a 

+  sec*0  B  (0  )T(0  )  exp{-2ik  [h  +  (1  +  c*A  / cos  0  )  ^*h  ]) 
aii  i  2  1  r  1  a  s 

-  sec*0  B  (0  )T(0  )  exp(-2ik  [h  +  (1  +  c*A  /cos*0  )X^*h  ]) 

X  }  a  2  r  j  1  2  i  2  2 


(46) 


where  Bx(9)  and  Ba(B)  are  defined  as  the  constant  density  Bayleigh 
reflection  coefficients  at  the  first  and  second  interfaces,  respectively. 
If  the  reference  velocity  is  the  velocity  in  the  first  medium,  then  A  =  0, 


and  the  density  artifact  is  given  by 


I  (-2k  )  =  - ; — : - - 
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+  sec*0  B  (0  )T(0  )  exp{-2ik  [h  +  (1  +  c*A  /cos*0  b  ]} 

S  S  1  1  r  x1  11  *  * 


sec  OjBj (OjJKOj)  exp{-2ik  [1^  +  (1  +  c1A1/cos  0^)  ^  hj) 


.  (47) 


The  expression  given  by  eqnstion  (47)  corresponds  to  the  spatial 
Fourier  transform  of  tbe  density  artifacts  plotted  in  Figures  25  through  32. 
To  understand  the  nature  of  these  artifacts,  ve  will  analyze  this  expression 
in  the  context  of  Born  theory.  Considering  the  first  term  in  this  equation, 
we  observe  that  this  expression  is  the  fixed  angle  analog  of  equation  (35). 
This  term  represents  the  artificial  density  variation  produced  by  the 
incomplete  reconstruction  of  the  bulk  modulus  at  the  first  interface. 
Reviewing  the  density  artifacts  plotted  in  the  aforementioned  figures,  we 
see  that  this  term  is  insignificant.  The  density  artifact  located  at  the 
first  interface  does  not  compromise  the  overall  accuracy  of  the 
reconstruction  since  the  input  angles  comply  with  the  Born  approximation. 
In  contrast,  the  second  and  third  terms  of  equation  (47)  combine  to  form  the 
deleterious  density  artifact  associated  with  the  second  interface. 

Before  examining  the  analytic  expression  for  the  second  interface 
density  artifact,  we  recall  some  assumptions  implicit  in  the  linear 
inversion.  The  primary  assumption  of  Born  inversion  is  that  individual 
scattering  interactions  do  not  significantly  affect  the  incident  wave  field. 
Thus,  the  total  wave  field  is  approximated  by  the  reference  wave  field,  and 
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the  inversion  models  all  wave  motion  as  the  wave  field  propagating  in  the 
reference  medium.  Consequently,  for  our  problem,  the  wave  field  propagating 
in  the  second  median  is  phase  shifted  with  respect  to  the  wave  field 
propagating  in  the  first  or  reference  median.  This  phase  shift  is  not 
accounted  for  by  the  inversion  and  contributes  to  inversion  errors.  In  onr 
analysis  of  single  interface  artifacts,  we  demonstrated  that  the  weak 
scatterer  assumption  is  satisfied  if  we  select  offset  data  which  minimize 
the  energy  of  the  reflected  wave  field.  As  ahown  above,  this  requirement  is 
not  sufficient  for  the  accurate  reconstruction  of  the  second  reflector.  To 
better  understand  inversion  errors,  we  must  consider  the  effect  of  the  [riiase 
shift  encountered  by  the  incident  field  at  the  first  interface. 

Returning  to  equation  (47),  we  observe  that  the  sum  of  the  second  and 
third  terms  of  this  expression  form  the  density  artifact  which  delineates 
the  second  interface.  In  general,  theae  two  terms  represent  two  delta 
functions  whose  amplitudes  are  similar  in  magnitude,  yet  opposite  in  sign. 
Furthermore,  these  delta  functions  are  characterized  by  different  phase 
shifts  which  affect  their  respective  locations  in  the  spatial  domain.  In 
fact,  these  two  phase  shifts  are  the  manifestation  of  the  phase  shifting  of 
the  incident  field  in  the  second  medium.  That  is,  if  A*  were  zero,  then  the 
wave  field  in  the  second  medium  would  be  in  phase  with  the  reference  wave 
field,  and  the  phase  shifts  which  characterize  these  two  delta  functions 
would  vanish.  As  a  result,  the  second  interface  artifact  would  be 
represented  by  a  single  delta  function  whose  amplitude  would  be  of  the  same 
order  as  that  of  the  single  interface  artifact.  Therefore,  as  the  velocity 
constrast  at  the  first  interface  decreases  in  accordance  with  the  Born 
approximation,  the  phase  shifts  inherent  in  equation  (47)  diminish,  and  the 
second  interface  artifact  decreases.  Upon  further  inspection,  we  conclude 


that  the  relative  amplitude  and  phase  differences  between  these  two  delta 
functions  are  dependent  on  the  inpnt  angles.  For  example,  as  6X  approaches 
6a»  both  the  phases  and  the  moduli  of  these  two  terms  become  equal,  and  the 
two  expressions  cancel.  Likewise,  as  the  values  of  8X  and  8a  diverge,  the 
two  delta  functions  become  more  dissimilar  in  amplitude  and  phase,  and  the 
artifact  increases.  Thus.  as  in  the  single  interface  situation,  the 
artifact  associated  with  the  second  interface  diminishes  as  6a  approaches 
®a.  This  second  interface  artifact  however  is  affected  in  both  amplitude 
■ad  phase.  Simply  stated,  it  is  the  phase  shifts  implicit  in  these 
expressions  which  are  responsible  for  the  deleterious  nature  of  the  second 
interface  artifacts. 

In  the  previous  discussion,  we  observed  that  since  the  incident  field 
is  phase  shifted  at  the  first  interface,  the  reconstructed  second  interface 
is  characterized  by  a  sum  of  two  phase  shifted  delta  functions.  These  phase 
shifts  are  represented  by 


where  cx  is  the  velocity  in  the  first  medium,  8j  is  the  incident  angle,  and 
Ax  is  the  velocity  contrast  at  the  first  interface.  As  mentioned  above, 
both  the  difference  between  these  phase  shifts  as  well  as  their  magnitudes 
influence  the  quality  of  the  reconstruction.  For  instance,  if  the  second 
term  of  equation  (48)  vanishes,  then  the  field  in  the  second  medium  is  in 
phase  with  the  incident  field,  and  inversion  results  are  improved.  Thus,  to 
rectify  inversion  errors,  we  should  minimize  the  size  of  this  term  in  the 
phase  shifts.  This  is  accomplished  by  considering  models  which  are 
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characterized  by  both  a  snail  velocity  contrast  at  the  first  interface,  and 


near  nornal  wave  vectors.  Considering  briefly  the  constant  density  models 
of  Fignres  25  throogh  27,  we  conclude  that  indeed  both  criteria  are 
necessary — only  the  nodel  of  Figore  27  where  0.1%  velocity  contrasts  are 
inverted  leads  to  an  insignificant  density  artifact  at  the  second  interface. 
Therefore,  by  adhering  to  the  weak  acatterer  assumption,  we  minimize  the 
phase  shift  of  the  incident  wave  field  which  enables  an  accurate 
reconstruction  of  the  second  reflector. 

In  examining  equation  (48),  we  note  that  only  the  strength  of  the  first 
velocity  contrast  affects  the  magnitudes  of  the  phase  shifts  implicit  in  the 
second  interface  density  artifact.  Consequently,  a  large  velocity  contrast 
at  the  first  interface  limits  the  accuracy  of  subsequently  reconstructed 
parameter  changes.  To  demonstrate  this  result,  we  consider  two  models  which 
consist  of  both  a  0.1%  and  a  10%  velocity  change.  In  these  two  models,  the 
order  of  the  velocity  changes  are  reversed  so  as  to  reveal  the  influence  of 
Ax  on  the  reconstructions  (see  Figures  30  and  31).  Turning  to  Figure  30 
first,  we  note  that  the  0.1%  velocity  change  at  the  first  interface  does  not 
seriously  impair  the  reconstruction  of  the  10%  velocity  change  at  the  second 
interface.  This  is  the  expected  result  as  is  trivial  in  this  case. 
Comparing  the  reconstructions  of  the  10%  velocity  changes  for  both  models, 
we  observe  that  the  corresponding  density  artifact  is  larger  in  Figore  30 
than  in  Figure  31.  Ve  also  note  that  the  actual  bulk  modulus  change  in 
Figure  30  is  slightly  antisymmetric.  (This  asymmetry  verifies  that  the 
actual  parameter  change  is  also  a  sum  of  two  band-limited  delta  functions.) 
Baaed  on  the  previous  discussion,  we  conclude  that  these  features  are 
accumuletive  manifestations  of  the  slight  phase  shift  incurred  at  the  first 
interface.  When  comparing  the  reconstructions  of  Figures  30  and  31 
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directly,  we  notice  that  they  are  not  simple  mirror  images  of  each  other. 
Since  A!  is  appreciable  for  the  model  of  Fignre  31,  the  reconstructed 
parameter  changes  associated  with  the  second  interface  are  noticeably 
incorrect  in  both  size  and  location.  Thus,  although  the  sizes  of  the 

reconstructed  parameter  changes  at  the  second  interface  are  weighted  by  the 
corresponding  reflection  coefficient,  the  overall  character  of  these 
variations  is  determined  by  the  phase  shifts  given  in  eqnation  (48). 

When  considering  the  general  features  of  double  interface 

reconstructions,  we  observed  that  errors  associated  with  the  Born 
reconstruction  of  the  first  interface  increase  with  depth.  These  errors 
affect  both  the  sizes  and  the  locations  of  the  reconstructed  parameter 
changes  delineating  the  second  interface.  By  examining  the  expression  for 
the  density  artifact  given  by  equation  (47),  we  can  understand  the  depth 
dependent  nature  of  the  error.  In  the  second  and  third  terms  of  equation 
(47),  we  notice  that  each  corresponding  phase  shift  is  weighted  by  the  depth 
of  the  second  reflector.  Ihos  as  h,  increases,  the  magnitudes  of  these 
phase  shifts  also  increase.  This  result  may  also  be  understood  in  the 

context  of  Born  theory.  In  Born  inverse  methods,  the  scattered  wave  field 
is  estimated  by  the  second  term  in  the  Born  series.  As  a  result,  the 

suitability  of  this  estimation  depends  upon  the  convergence  of  the  higher 
order  terms.  If  this  series  converges  slowly,  then  the  first  Born 
epproximation  inadequately  describes  the  scattered  wave  field.  In  any  case, 
since  the  higher  order  terms  in  the  Born  series  are  integral  expressions, 
the  error  associated  with  their  convergence  increases  as  the  limits  of 
integretion  increase.  That  is,  errors  implicit  in  the  lineer  inversion 
accumulate  as  we  integrate  the  wave  field  over  greater  and  greater  depths. 
In  conclusion,  we  maintain  that  the  predicted  sizes  and  locations  of  the 
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second  parameter  changes  are  determined  by  both  the  phase  shifting  of  the 
incident  wave  and  the  depth  of  the  second  interface. 

Inversion  Artifacts  and  the  Reference  Velocity 

In  the  previoos  reconstructions,  the  reference  velocity  is  chosen  to  be 
the  velocity  in  the  first  medium.  As  a  result,  the  wave  field  in  the  second 
median  is  phase  shifted  with  respect  to  the  wave  field  in  the  reference 

medium,  and  this  causes  inversion  errors  to  increase  with  depth.  In  Born 
inversion  however,  we  need  not  choose  the  reference  velocity  to  be  the 
velocity  in  the  first  median.  The  specific  choice  of  s  reference  velocity 
depends  on  the  goals  of  the  inversion.  In  this  section,  we  examine  the 
influence  of  the  reference  velocity  on  inversion  artifacts.  To  do  this,  we 
invert  a  constant  density  model  for  different  reference  velocities  and 
interpret  the  results  in  terms  of  equation  (46).  As  mentioned  above, 

equation  (46)  represents  the  double  Interface  density  artifact  for  an 

arbitrary  reference  velocity.  From  this  expression,  we  observe  that  the 

reconstruction  of  each  interface  is  the  summation  of  two  phase  shifted  delta 
functions.  Thus,  if  the  reference  velocity  is  neither  ca  nor  ca.  then  each 
reconstructed  interface  is  characterized  by  the  accumulative  errors  produced 
by  these  phase  shifts.  As  a  result,  not  only  are  both  the  sizes  and 
locations  of  the  parameter  changes  incorrect  for  each  interface,  but  also 
the  extent  of  these  errors  are  determined  by  the  depths  of  the  interfaces. 
As  before,  we  note  that  the  strengths  of  these  phase  shifts  depend  on  the 
angle  of  incidence,  the  velocity  contrast,  and  the  reference  velocity.  We 
also  recall  that  these  phase  shifts  which  spatially  separate  the  delta 
functions  only  exist  if  the  wave  field  above  the  reconstructed  interface  is 
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phase  shifted  with  respect  to  the  reference  wave  field. 


In  Figure  32,  we  present  the  reconstruction  of  a  constant  density  model 
with  the  reference  velocity  equal  to  the  velocity  in  the  first  medium.  In 
this  example,  the  reconstruction  of  the  first  interface  is  essentially 
correct  while  the  reconstruction  of  the  second  interface  suffers  from 
accumulative  phase  errors.  Considering  the  first  reconstructed  reflector, 
we  observe  that  the  locations  of  the  parameter  changes  are  correct/  Born 
theory  however  predicts  slight  errors  in  the  reconstruction  of  their 
corresponding  amplitudes.  For  our  purposes,  we  view  this  reconstruction  as 
the  prototype  for  both  interfaces  in  all  subsequent  examples.  To  provide  a 
proper  comparison  between  this  reconstruction  and  other  relevant 
reconstructions,  we  list  the  amplitudes  of  these  variations  in  Table  1  (see 
page  55)  . 

In  Figure  33,  we  reconstruct  the  model  of  Figure  32  with  the  reference 
velocity  equal  to  the  layer  velocity.  Consequently,  the  wave  field  in  the 
layer  is  in  phase  with  the  reference  wave  field  since  Ax  =  0.  In  this 
example  however,  the  wave  field  in  the  first  medium  is  phase  shifted  with 
respect  to  the  reference  wave  field  (A  h  0) .  We  conclude  that  the  artifact 
associated  with  the  first  interface  is  characterized  by  two  phase  shifted 
delta  functions  while  that  of  the  second  interface  consists  of  only  one 
delta  function.  Thus,  as  the  wave  field  propagates  out  of  phase  to  the 
first  reflector,  the  error  inherent  in  the  Born  reconstruction  of  this 
reflector  grows  in  size.  From  Figure  33,  we  observe  that  this  accumulative 
error  affects  both  the  sizes  and  the  locations  of  the  parameter  changes. 
Comparing  the  parameter  amplitudes  of  this  reconstruction  with  their 
counterparts  in  Figure  32  (see  Table  1),  we  conclude  that  both  parameter 
changes  have 
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increased  in  size.  In  fact,  there  is  now  a  noticeable  density  artifact 
located  at  the  first  interface.  Considering  the  second  reflector,  we 
observe  that  the  error  in  this  reconstruction  is  similar  to  that  of  the 
first  interface  (see  Table  1). 


FIGURE 

E  max 

X 

p  max 

K  i 

E  max 

2 

P  j®ax 

32 

0.093 

— 

0.140 

0.108 

33 

0.103 

0.005 

0.103 

0.005 

35 

0.095 

0.002 

0.095 

0.002 

Table  1.  List  of  the  maximum  parameter  amplitudes  for  Figures  32,  33, 
and  35. 

Since  the  wave  field  in  the  layer  it  in  phase  with  the  reference  wave  field, 
the  errors  established  at  the  first  interface  propagate  to  the  second 
reflector  bat  do  not  grow  in  size.  Thus,  the  location  error  at  the  second 
interface  is  identical  to  that  at  the  first  interface.  There  are  however 
slight  amplitade  errors  associated  with  the  incomplete  reconstruction  of  the 
second  reflector.  To  further  illustrate  that  these  inversion  errors  do  not 
increase  within  the  layer,  we  consider  this  model  for  a  deeper  second 
reflector  (see  Figure  34).  In  comparing  Figures  33  and  34,  we  conclude  that 
the  inversion  errors  associated  with  the  second  interface  are  identical  to 
those  of  the  first.  Therefore,  by  establishing  the  reference  velocity  as 
the  layer  velocity,  we  have  distributed  the  inversion  error  between  both 

reflectors,  and  now  these  errors  are  sensitive  to  the  depth  of  the  first 

refl actor . 

As  mentioned  above,  if  the  reference  velocity  is  neither  c,  nor  c,, 

then  both  reflectors  are  reconstructed  by  a  sum  of  two  phase  shifted  delta 

functions.  As  the  depths  of  these  reflectors  increase,  the  corresponding 


phase  discrepancies  also  increase,  and  the  inversion  errors  grow  in  size. 
These  observations  follow  directly  from  equation  (46).  To  demonstrate  this 
result,  we  reconstruct  the  previous  constant  density  model  for  a  reference 
velocity  midway  between  cx  and  c,;  in  Figure  35,  we  present  this  example 
with  cr  =  5025  ft/sec.  Comparing  this  reconstruction  with  that  of  Figure  33, 
we  observe  that  the  inversion  errors  associated  with  the  first  interface 
have  diminished.  By  reducing  the  reference  velocity,  we  have  reduced  phase 
shifts  in  the  first  two  terms  of  equation  (46),  thereby  reducing  the 
accumulative  errors.  Thus,  both  the  density  artifact  and  the  location  error 
decrease  at  the  first  interface.  From  Table  1,  we  also  note  that  the 
amplitude  of  the  bulk  modulus  variation  is  also  reduced.  In  fact,  it  is 
close  to  its  counterpart  in  Figure  32.  Hnlike  the  model  of  Figure  33,  the 
wave  field  in  the  layer  of  Figure  35  is  phase  shifted  with  respect  to  the 
reference  wave  field.  As  a  result,  the  inversion  errors  produced  at  the 
first  interface  not  only  propagate  with  depth,  but  also  they  increase  with 
depth.  To  demonstrate  this  feature,  we  invert  the  model  of  Figure  35  for  a 
deeper  second  reflector  (see  Figure  36).  Clearly,  the  reconstruction  of  the 
second  reflector  in  Figure  36  is  characterized  by  greater  inversion  errors-- 
the  predicted  parameter  variations  are  both  larger  and  shallower.  In 
conclusion,  we 

observe  that  this  intermediate  reference  velocity  effectively  distributes 
the  inversion  error.  In  fact,  this  reconstruction  is  superior  to  that  of 
Figure  33  as  the  locations  and  the  amplitudes  of  both  velocity  changes  are 
more  accurate. 

In  the  above  examples,  we  demonstrate  that  by  adjusting  the  reference 
velocity,  we  may  distribute  the  inversion  errors  between  the  reconstruction 
of  each  interface.  Depending  on  the  goals  of  the  inversion,  this  procedure 
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might  be  advantageous.  For  example,  when  comparing  Figures  32  and  35,  we 
might  conclude  that  Figure  35  is  a  more  accurate  reconstruction  of  the  input 
model.  Even  though  all  amplitudes  and  locations  are  incorrect,  this 
reconstruction  is  not  dominated  by  the  error  associated  with  the  second 
reflector.  We  might  however  wish  to  reconstruct  the  first  interface  with 
precision.  In  that  case,  the  reconstruction  of  Figure  32  woold  be  superior 
to  that  of  Figure  35.  Therefore,  since  errors  are  an  inevitable  consequence 
of  constant  background  inversion,  we  choose  a  reference  velocity  so  as  to 
minimize  the  errors,  deemed  most  deleterious  to  our  goals. 
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CONCLUSIONS 


In  this  thesis,  we  implement  the  two-parameter,  constant  background. 
Born  inversion  method  of  Clayton  and  Stolt  (1981)  for  acoustic  synthetic 
data  which  represent  single  and  double  interface  models.  Since  these  data 
are  exact,  we  are  able  to  compare  reconstructed  acoustic  profiles  with  their 
analytic  expressions.  We  direct  this  analysis  toward  understanding  the 
error  implicit  in  the  approximate  inversion.  As  a  result,  we  devise  an 
implementation  which  provides  more  reliable  interpretations  of  Born 
reconstructions. 

From  our  analysis,  we  observe  that  inversion  errors  are  minimized  if 
the  offset  data  is  preconditioned  in  accordance  with  the  Born  approximation. 
Specifically.  we  reduce  the  energy  of  the  reflected  wave  field  by 
eliminating  the  offset  data  which  are  characterized  by  large  angles  of 
incidence.  These  deleterious  data  however  provide  the  greatest  variation  in 
the  reflection  coefficient.  As  pointed  out  by  Clayton  and  Sto.t  (1981),  it 
is  the  angular  dependence  of  the  reflection  coefficient  which  allows  the 
separation  of  the  components  of  reflectivity.  Therefore,  inherent  in  our 
method  is  a  contradiction  of  objectives:  the  linear  inversion  requires  the 
elimination  of  data  which  enhance  the  recovery  of  both  acouttic  parameters. 

By  selecting  a  fixed  angle  implementation,  we  are  able  to  separate 
bandwidth  considerations  from  data  preconditioning  requirements. 
Consequently,  the  flexibility  afforded  by  this  implementation  leads  to 
stable  results.  In  fact,  for  single  interface  models,  the  reconstructed 
profiles  are  almost  exact.  For  double  interface  models  however,  the 
reconstructions  suffer  from  sizeable  Born-related  artifacts.  To  understand 


the  natore  of  these  artifacts,  we  examine  their  analytic  expressions.  Since 


the  reconstructed  density  variation  is  less  reliable  than  its  corresponding 
bulk  modulus  variation,  we  focus  our  analysis  on  double  interface  density 
artifacts.  As  a  result,  we  test  the  inversion  for  the  input  models  which 
provide  the  least  favorable  results. 

In  our  analysis  of  double  interface  density  artifacts,  we  observe  that 
the  reference  velocity  influences  the  interpretabil i ty  of  the 
reconstructions.  In  general,  the  reference  velocity  affects  the  phase 
shifts  inherent  in  the  modeled  wave  field.  As  evidenced  by  our  examples, 
this  phase  shifting  causes  inversion  artifacts  to  increase  with  depth  —  these 
accumulative  phase  errors  seriously  can  compromise  the  usefulness  of  this 
inversion  technique.  We  therefore  select  a  reference  velocity  which  reduces 
the  errors,  considered  most  harmful  to  our  inversion  objective. 

This  investigation  establishes  that  this  two-parameter,  constant 
background,  linear  inversion  cannot  accurately  reconstruct  two  consecutive, 
real  earth,  parameter  changes,  even  under  idealized  conditions.  The  data 
requirements  preclude  the  applicability  of  the  Born  approximation  for 
recovering  the  acoustic  parameter  changes,  typical  of  geophysical  problems. 
Ihus,  this  inversion  accurately  reconstructs  only  weak  scattering 
interactions.  To  improve  the  usefulness  of  this  technique,  we  could 
implement  the  inversion  for  a  variable  reference  velocity.  If  accurate 
migration  velocities  are  obtained,  this  approach  will  reduce  the  relative 
sizes  of  parameter  changes,  and  hence  validate  the  weak  scatterer 
assumption.  As  a  result,  accumultive  phase  errors  will  be  reduced. 

Applying  the  refinement  algorithm  developed  by  Hagin  and  Cohen  (1984) 
to  this  two-parameter  procedure  would  reduce  the  accumulative  errors, 
implicit  in  the  linear  inversion.  According  to  this  refinement  algorithm. 
Born-related  errors  are  recursively  eliminated  from  constant  background 
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reconstructions. 


Consequently,  this  procedure  provides  an  improved  estimate 


of  snbsorface  parameter  variations.  This  method  however  requires  an 
accurate  estimate  of  real  earth  parameter  changes  from  measurements  of  their 
reconstructed  counterparts.  To  insure  accurate  estimates,  a  method  for 
separating  artifacts  from  true  parameter  variations  would  have  to  be 
developed.  Perhaps  however,  the  single  interface  artifact  identification 
procedure  presented  by  Weglein,  Violette,  and  Keho  (1985)  could  be  used  to 
test  each  successive  reflector  for  inversion  artifacts.  Since  the 
aforementioned  refinement  scheme  would  eliminate  accumulative  errors,  this 
identification  procedure  could  be  applied  with  confidence  at  each  interface. 
We  note  that  this  method  would  be  most  useful  for  noisy  synthetic  or  field 
data  sets.  Implementing  the  refinement  procedure  in  conjunction  with  this 
artifact  identification  method  would  be  computationally  more  efficient  than 
implementing  a  variable  background  inveraion.  Furthermore,  these  methods 
require  only  an  accurate  estimate  of  the  velocity  in  the  first  medium, 
whereas  variable  background  inversion  requires  an  accurate  estimate  of  all 
migration  velocities.  Therefore,  this  postprocessing  refinement  algorithm 
might  improve  two-parameter  reconstructions  by  efficiently  reducing  the 
effects  of  accumulative  errors. 
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APPENDIX  A:  TKANSFOHED  G KEEN’S  FUNCTIONS 


In  this  Appendix,  we  present  the  Fourier  transforms  of  the  analytic 
expressions  for  the  Green's  functions  in  equation  (10).  These  Green's 
functions  are  the  resolven'  for  the  constant  background,  two-dimensional 
wave  operator.  If  pr  and  are  assumed  constant,  Lr  becomes 


L  = 
r 


V  + 


( A-l ) 


The  Green's  function  representing  wave  propagation  from  the  scatterer 
location  to  the  geophone  location  solves 


-p  8(x 
r  S 


x')6(z’) 


( A-2) 


Introducing  the  spatial  Fourier  transform  over  the  receiver  coordinates  as 


6  (k  ,  0  |x  ’,  z '/«) 
r  xg  1 


G  (x  ,  0  lx',  z  ’/«) 
r  g  1 


exp("ikxgxg)  ' 


( A-3 ) 


we  apply  this  transform  to  equation  (A-2)  to  obtain  the  following  ordinary 
differential  equation  for  the  function  &r: 
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I  — -  +  k  G  (k  ,0  x',z'/w)  =  - —  exp(-ik  x' 

l  dz1  z»  J  r  xg  fa  xg 


)  6  ( z  ' ) 


(  A  —  4 ) 


where  k 


iL  _  k 
■:  ** 


Two  linearly  independent  solntions  to  the  homogeneous  form  of  ( A— 4)  are 
given  by 


G  =  - — 1  exp(+ik  |z  -  z'l)  exp(-ik  x  ) 

r  fa  z«'  '  x«  « 


( A-5) 


where  represents  the  outgoing  eolation  and  Gf-  represents  the  incoming 

solution.  For  our  purposes,  we  are  interested  only  in  the  outgoing  solution 
and  have  only  to  determine  the  quantity  A+. 

Ve  determine  A4  by  imposing  two  conditions  at  the  scatterer  depth  which 
characterize  the  nature  of  the  scatterer.  These  conditions  are:  (1)  that 
the  function  be  continuous,  and  (2)  that  its  first  derivative  have  a  jump 
equal  to  the  source  term  over  an  interval  containing  z'.  These  conditions 
yield: 


+  i 

A  -  yy—  / 

zg 


( A— 6 ) 


and.  at  the  receiver  location,  the  Green's  function  is  given  by 


<  A  -  7 ) 


&  a  »o|x’,z'/w) 
r  xg  1 


lpr 

-j —  exp[-i(k  x' 
1 2  n  8 


|z  '  | ) ] / 2kz^ 


Ve  now  perform  a  similar  analysis  to  determine  the  transformed  Green’s 
fnnction  which  characterizes  propagation  from  the  soorce  location  to  the 
scatterer  location.  The  only  difference  between  this  derivation  and  the 
previons  one  is  that  the  Fourier  transform  over  the  soorce  coordinates  has 
the  opposite  sense  as  equation  (A-3).  This  difference  results  in  a  sign 
change  on  the  argument  of  the  exponential: 


ff+(x\z'|k  ,  Ofoi) 
r  1  xs 


ipr 

—  exp[-i(kxfx'  +  ktJz'|))/2kzs  . 
12n 


where 


k 

zs 


( A-8) 


Ve  now  make  the  observation  that  the  parameter  variations  are  zero  for  z 
less  than  zero.  This  result  allows  us  to  drop  the  absolute  value  signs  in 
equations  (A-7)  and  (A-8)  before  substitution  into  equation  (10) 
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APPENDIX  B:  ONE-DIMENSIONAL  FOE*  AID  EQUATION 


In  this  Appendix,  we  derive  the  one-dimensional  forward  equation 
from  the  two-dimensional  forward  eqnation.  This  is  accomplished  by 
assnming  that  the  probed  medium  is  independent  of  both  transverse 
directions.  As  a  result,  the  observed  offset  data  is  only  a  function 
of  xg  -  xg.  Explicitly  writing  the  spatial  Fourier  transform  of  the 
left  hand  side  of  equation  (13) 


x  ,o>)exp[-i(k  x  -  k  x  )],  ( B— 1 ) 

s  xg  g  xs  s 


we  express  this  result  in  terms  of  midpoint,  half-offset  coordinates: 


jl  =  (x  -  x  )/2  ,  x  «  ( x  +  x  )/2  (B-2) 

h  g  s  mgs 


to  obtain 


D(k  ,k  ,  u>) 

xg  xs 
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2n 


2dx^  D(2x^,a>) 


•  exp[-ik  (x.  +  x  )]  exp[ik  (x  -  x_  )  ] 

xg  n  m  xs  m  n 


(B-3) 
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where  dxgdxs  =  2dxj,dxni 
Rewriting  equation  (B-3) 


D(k 


CD 

’kxs'w)  =  57  { 


2dx  D(  2x  ,  (o)  exp[-i  (k  +  k  )x .] 
h  h  xg  xs  h 


dx  exp[-i(k  -  k  )x  ]  , 

m  xg  xs  m 


( B—  4 ) 


we  recognize  both  integrals  as  Fourier  transforms: 


B(k  ,k  ,(o)  =  6(k  -  k  ) B[ ( k  -  k  ) / 2 ,  w] 

xg  xs  xg  xs  xg  xs 


(B-5) 


Turning  to  the  right  band  side  of  equation  (13).  we  write  the 
transformed  model  variations  as  Fourier  transforms 


•  exp[-i(k  -  k  ) ] x  exp[i(k  +  k  )z]  (B-6) 

xg  xs  xs  zg 

where  aj  represents  both  parameter  variations.  For  a^  independent  of 
x.  this  result  becomes 


B-7) 


*.  ( k  -  k  ,-k  -  k  )  =  6(k  -  k  )  I .  ( -k  -  k  )  . 

1  Xg  XS  Zg  ZS  Xg  XS  1  zg  zs 


Thus,  for  a  one-dimensional  median,  the  linear  eqaations  are: 
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Integrating  both  sides  of  equation  (B-8)  with  respect  to  kxs»  »nd 
observing  that  kxg  =  kxs  and  kxg  =  kx#,  we  find 


i  ,«)  =  — T-  f  k*I  (-2k  )  +  (k*  -  k* )  a  (-2k  )  1 

z  ,vl  l  r  »  z  z  x*  zj 


S ( w)  .  (B-9) 


By  the  dispersion  relation,  the  right  hand  side  of  (B-9)  is  expressed 
only  as  a  function  of  k  yielding  equation  (14). 
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APPENDIX  C:  LEAST-SQUARES  SOLUTION  OVER  TEMPORAL  FRBQUENC! 


The  forward  equation  expressed  by  equation  (15)  can  be  solved  for  the 
mediant  parmeters  by  applying  1  ea  s  t -sq  uar  e  s  principles.  In  general,  the 
least  sqoares  criterion  maintains  that  we  minimize  the  square  of  the 
difference  between  the  desired  and  the  actual  expressions.  In  our  case,  the 
desired  earth  response  is  the  data  itself.  The  actual  response  however  is 
the  right  hand  side  of  equation  (15).  Thus,  the  least  sqoares  error  is: 


(C-l) 


where  the  summation  is  taken  over  a  specified  range  of  frequency  values. 
This  error  is  minimized  by  setting  its  partial  derivatives  with  respect  to 
each  parameter  equal  to  zero: 


I7  =  2}  [  D.  -  A  .a  -  AT  lAm|  =  0 
f  l  j  »j  *  J  ^ 

J  ( C-2  > 

4—  =2^[d  -At  -  A  f  I A  =  0 

^  {  l  j  »j  *  *  J  »J 


The  resulting  system  of  equation  is  known  as  the  normal  equations  given  by 
equations  (18).  The  solution  of  the  normal  equations  yields  the  least- 
squares  determination  of  the  medium  parameters. 


APPEND!!  D:  SINGLE  INTERFACE  LINE  SOURCE  ACOUSTIC  DATA 


In  this  Appendix,  we  reduce  the  point  source  expression  for  wave 
propagation  in  two  semi-infinite  media  to  the  line  source  result.  On 
page  96  of  Ewing,  Jardetzky,  and  Press  (1957),  the  scattered  wave  field 
in  the  first  medium  is  given  by 


f  6*x  -  V, 

Dg (  r,  go)  =  I  dk  ^ — T~y~  ^exPl“i|,1^z  +  h)]J#(kr)k}/^J  (D-l) 


where  6  =  px/pj,  an^ 

interface.  In  this 
coordinates  and  thus: 


h  is  the  distance  of 
expression,  r  and  z 


the  source  above  the 
represent  cylindrical 


= 


(  D  -  2  ) 


where  k 


For  layered  media,  the  above  point  source  result  has  circular 
symmetry,  and  we  recognize  equation  (D-l)  as  the  Foorier-Bessel 
transform.  Rewriting  equation  (D-l)  as  t  Fourier  transform  in 


rectangular  coordinates,  the  field  at  z  =  h  becomes 


D  (x,  y,  u,)  = 
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j  1 


dk 


sy,  -  i/2 


y  f>v2  + 

(exp(-2^ih)]/*»1  .  exp[-i(k^x  +  k^y)] 


( D-3 ) 


Since  the  probed  medium  is  assumed  independent  of  the  y-direction,  we 
may  integrate  over  this  direction  to  obtain  the  line  source  result 


D  (  x,  oj)  =  2  n 


£.  y  -  y 

dkx  bV~~+  V  ^exP(-2*,1h,/*'J  exp(-ikxx) 
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where  now 


"x  =  1 
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Fourier  transforming  equation  (0-4),  we  find 


D  ( k  ,  <d)  = 


6 1 ,  - 


(D-6) 
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For  interpretational  purposes,  Cohen  and  Bleistein  (1982)  suggest 
that  it  is  more  desirable  to  process  data  for  the  normal  derivative  of 
the  variation  rather  than  the  variation  itself.  For  our  purposes,  this 
derivative  is  achieved  by  multiplying  equation  (D-6)  by  yx  which  leaves 
equation  (19) . 


APPENDIX  E:  DERIVATION  OF  THE  DOUBLE  INTERFACE  FOKVARD  MODEL 


The  purpose  of  this  Appendix  is  to  derive  the  expression  for  the 
scattered  wavefield  generated  by  the  two-interface  acoustic  model.  The 
specifics  of  this  model  are  illustrated  in  Figure  3.  The  method  of 
determination  is  to  solve  the  system  of  equations  which  satisfy  the  boundary 
conditions  at  each  interface  for  the  displacement  potential  representing  the 
upgoing  wave  field  in  the  first  medium.  The  details  of  the  procedure  are 
similar  to  examples  presented  by  Ewing,  Jardetzky,  and  Press  (1957). 

Omitting  the  time  factor,  we  can  express  the  potential  in  the  first 
medium  as 


03 

‘A 


J  (kr)  exp(-yiz)  + 


CO 
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J#(kr)  expf-^z) 


(E-l) 


where  Qx  is  the  quantity  to  be  determined.  All  other  quantities  are  defined 
in  Chapter  2  as  well  as  in  Ewing,  Jardetzky,  and  Press.  The  first  term  in 
equation  (E-l)  represents  the  direct  wave  while  the  second  term  corresponds 
to  the  upgoing  or  scattered  wave  field.  The  potential  in  the  second  medium 
is 


'■-R 


J#(kr)  exp(-yjz)  + 


CD  * 
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J#(kr)  exp(yjZ) 


( E-2 ) 


These  two  terms  correspond  to  the  summation  of  the  upward  and  the  downward 
traveling  waves  which  reflect  many  times  within  the  layer.  Finally,  the 
potential  in  the  third  medium  is  given  by 


K 


dk  —  Ja(kr) 


exp(yf  z) 


(E-3) 


which  represents  the  downgoing  wsve  field  in  the  lower  half-space. 

The  potentials  and  0a  most  satisfy  the  boundary  conditions: 

P.P.  =  P.P.  »  77-  =  37-  (E-4) 

XX  22  0  Z  OZ 

at  the  first  interface.  These  conditions  express  the  notion  that  both  the 
pressure  and  the  normal  displacement  are  continuous  across  the  interface. 
Similarly,  the  potentials  in  the  second  and  third  media  must  satisfy  the 
boundary  conditions: 


dp  dp 

P2P2  =  P,P,  .  =  ■57-  (E_5) 

at  the  second  interface.  These  conditions  lead  to  a  system  of  four  linear 
equations  expressed  by 


where 
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As  before,  we  are  interested  in  the  plane-wave  reflection  coefficient, 
and  this  quantity  corresponds  directly  to  Q14  Solving  for  Q1#  we  obtain 
equation  (21).  An  expression  similar  to  equation  (21)  is  presented  by  Frisk 
(1979)/  in  his  paper.  Frisk  discusses  many  characteristics  of  the  plane- 
wave  reflection  coefficient  for  single  and  double  interface  models. 

By  representing  the  denominator  of  equation  (21)  as  an  infinite  sum,  we 
are  able  to  decompose  the  forward  model  into  its  component  reflected  events. 
Thns,  equation  (21)  becomes 


CD 

Q”  =  exp(-2yih1)  |rx  +  8Jexp(-2yjhJ)  ]  )  ( -1 )  “(R  ,R, )  Bexp(-2n*rih ,) .  (E-7 
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Carrying  ont  a  portion  of  the  summation,  and  collecting  terms,  we  obtain 


Qx  =  R1exp(-2*/1h1)  +  R  ,  ( 1  -  R*)exp[-2(jraha  +  |fah#>] 

+  RaR*(l  -  Ra)ezp[-2(yaha  +  2*/JhJ>]  +  ...  <  E-8) 


where:  the  first  term  corresponds  to  the  reflection  from  the  first 
interface,  the  second  term  corresponds  to  the  reflection  from  the  second 
interface,  and  the  third  term  corresponds  to  the  first  multiple  reflection 
from  within  the  layer.  By  writing  the  data  in  this  form,  we  are  able  to 
establish  a  forward  model  which  excludes  multiple  reflections.  Not  only 
does  such  an  expression  provide  more  favorable  data  for  the  linear 
inversion,  but  also  it  yields  a  less  complicated  analytic  expression.  This 
result  is  given  by  equation  (25). 


Figure  1.  Diagram  of  the  tingle  interface  experiment. 


INPUT  MODEL 


Cl  = 

500D. 0 

D1  =  1.0000 

C2  = 

5500.0 

D2  =  1 . 0000 

FI  = 

60.0 

F2  =  62.0 

REFLECTION  COEFF. 


DENSITY  VARIATION  max  =  l.n 


BULK  MODULUS  VARIATION  imrx  =  1.3) 


1 

1 

—  v  r*\  /A  rA  r 

- 1 - - — 1 - 

v  y  U 

_ 1 _ 

pvv^- — - 

_ 1 _ 1 _ 

0  see  100a  1500  2000  ssez 

DEPTH  IFEET) 


Figure  5.  Fixed  frequency,  incretting  velocity,  constmt  demity 
•ingle  interface  reconstruction  for  ell  pre-critical  data. 
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Figaro  9.  Fixed  fraqotney,  iacreasiag  velocity.  contttat  iapodance 
eiagle  interface  reconstruction  for  all  pro-critical  data. 
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Figure  11.  Fixed  frequency,  decreasing  velocity,  conetent  density, 
single  interface  reconstrnction  for  broadband  date. 
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Figure  18.  Fixed  freqoency,  constant  density,  doable  interfeee  reeon- 
strnction  for  0.1%  Telocity  changes. 
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20.  Fixed  frequency,  conetent  bulk  ■odnlnt.  doable  interface  re 
construction  for  1%  velocity  changes. 
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21.  Fixed  frequency,  conetent  balk  Bodnias,  doable  interface 
construction  for  0.11  Telocity  chanfee. 
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Flgara  28.  lacoastractloa  of  tha  aodal  of  Figera  26  bat  for  a  shallower 
aacoad  raflactor. 
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Figaro  SO.  Fiiad  uglii  eoaitut  dmttj,  doable  Interface,  rocoo 
atractloD  for  0.11  and  101  ailotitj  chaagaa. 
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Figara  36.  laconatrectioa  of  Figara  35  bat  for  a  daeper  aaeood  ra 
flaetor . 
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